NASA/TM— 20 10-21 6766 


AIAA-20 10-543 



A New Approach for Constructing 
Highly Stable High Order CESE Schemes 


Sin-Chung Chang 

Glenn Research Center, Cleveland, Ohio 


August 2010 


NASA STI Program ... in Profile 


Since its founding, NASA has been dedicated to the 
advancement of aeronautics and space science. The 
NASA Scientific and Technical Information (STI) 
program plays a key part in helping NASA maintain 
this important role. 

The NASA STI Program operates under the auspices 
of the Agency Chief Information Officer. It collects, 
organizes, provides for archiving, and disseminates 
NASA’s STI. The NASA STI program provides access 
to the NASA Aeronautics and Space Database and 
its public interface, the NASA Technical Reports 
Server, thus providing one of the largest collections 
of aeronautical and space science STI in the world. 
Results are published in both non-NASA channels 
and by NASA in the NASA STI Report Series, which 
includes the following report types: 

• TECHNICAL PUBLICATION. Reports of 
completed research or a major significant phase 
of research that present the results of NASA 
programs and include extensive data or theoretical 
analysis. Includes compilations of significant 
scientific and technical data and information 
deemed to be of continuing reference value. 
NASA counterpart of peer-reviewed formal 
professional papers but has less stringent 
limitations on manuscript length and extent of 
graphic presentations. 

• TECHNICAL MEMORANDUM. Scientific 
and technical findings that are preliminary or 
of specialized interest, e.g., quick release 
reports, working papers, and bibliographies that 
contain minimal annotation. Does not contain 
extensive analysis. 

• CONTRACTOR REPORT Scientific and 
technical findings by NASA-sponsored 
contractors and grantees. 


• CONFERENCE PUBLICATION. Collected 
papers from scientific and technical 
conferences, symposia, seminars, or other 
meetings sponsored or cosponsored by NASA. 

• SPECIAL PUBLICATION. Scientific, 
technical, or historical information from 
NASA programs, projects, and missions, often 
concerned with subjects having substantial 
public interest. 

• TECHNICAL TRANSLATION. English- 
language translations of foreign scientific and 
technical material pertinent to NASA’s mission. 

Specialized services also include creating custom 

thesauri, building customized databases, organizing 

and publishing research results. 

For more information about the NASA STI 

program, see the following: 

• Access the NASA STI program home page at 
http://www.sti.nasa.gov 

• E-mail your question via the Internet to help@ 
sti.nasa.gov 

• Fax your question to the NASA STI Help Desk 
at 443-757-5803 

• Telephone the NASA STI Help Desk at 
443-757-5802 

• Write to: 

NASA Center for AeroSpace Information (CASI) 
7115 Standard Drive 
Hanover, MD 21076-1320 


NASA/TM— 20 10-21 6766 


AIAA-20 10-543 



A New Approach for Constructing 
Highly Stable High Order CESE Schemes 


Sin-Chung Chang 

Glenn Research Center, Cleveland, Ohio 


Prepared for the 

48th Aerospace Sciences Meeting 

sponsored by the American Institute of Aeronautics and Astronautics 
Orlando, Florida, January 4-7, 2010 


National Aeronautics and 
Space Administration 


Glenn Research Center 
Cleveland, Ohio 44135 


August 2010 


This work was sponsored by the Fundamental Aeronautics Program 
at the NASA Glenn Research Center. 


Level of Review. This material has been technically reviewed by technical management. 


Available from 


NASA Center for Aerospace Information 
7115 Standard Drive 
Hanover, MD 21076-1320 


National Technical Information Service 
5301 Shawnee Road 
Alexandria, VA22312 


Available electronically at http://gltrs.grc.nasa.gov 


A NEW APPROACH FOR CONSTRUCTING HIGHLY STABLE HIGH ORDER CESE SCHEMES 


Sin-Chung Chang 

NASA Glenn Research Center, Cleveland, OH 44135 
E-mail: sin-chung.chang@nasa.gov 

Abstract 

A new approach is devised to construct high order CESE schemes which would avoid the common 
shortcomings of traditional high order schemes including: (a) susceptibility to computational instabilities; 
(b) computational inefficiency due to their local implicit nature (i.e., at each mesh points, need to solve a 
system of linear/nonlinear equations involving all the mesh variables associated with this mesh point); (c) 
use of large and elaborate stencils which complicates boundary treatments and also makes efficient parallel 
computing much harder; (d) difficulties in applications involving complex geometries; and (e) use of problem- 
specific techniques which are needed to overcome stability problems but often cause undesirable side effects. 
In fact it will be shown that, with the aid of a conceptual leap, one can build from a given 2nd-order CESE 
scheme its 4th-, 6th-, 8th-,... order versions which have the same stencil and same stability conditions 
of the 2nd-order scheme, and also retain all other advantages of the latter scheme. A sketch of 
multidimensional extensions will also be provided. 

1. Introduction 

The space-time conservation element and solution element (CESE) method is a high-resolution and 
genuinely multidimensional method for solving conservation laws [1-63]. Its nontraditional features include: 
(i) a unified treatment of space and time; (ii) the introduction of conservation elements (CEs) and solution 
elements (SEs) as the vehicles for enforcing space-time flux conservation; (iii) a time marching strategy 
that has a space-time staggered stencil at its core and, as such, fluxes at an interface can be evaluated 
without using any interpolation or extrapolation procedure (which, in turn, leads to the method’s ability 
to capture shocks without using Riemann solvers); (iv) the requirement that each scheme be built from a 
nondissipative core scheme and, as a result, the numerical dissipation can be controlled effectively; and (v) 
the fact that mesh values of the physical dependent variables and their spatial derivatives are considered as 
independent mesh variables to be solve for simultaneously. Note that CEs are non-overlapping space-time 
subdomains introduced such that (i) the computational domain can be filled by these subdomains; and (ii) 
flux conservation can be enforced over each of them and also over the union of any combination of them. On 
the other hand, SEs are space-time subdomains introduced such that (i) the boundary of each CE can be 
divided into component parts with each of them belonging to a unique SE; and (ii) within a SE, any physical 
flux vector is approximated using simple smooth functions. In general, a CE does not coincide with a SE. 

Without using flux-splitting, dimensional-splitting, mesh-alignment or other special techniques, since its 
inception in 1991 [1], the unstructured-mesh compatible 2nd-orcler CESE method has been used to obtain 
numerous accurate ID, 2D and 3D steady and unsteady flow solutions with Mach numbers ranging from 
0.0028 to 10 [15]. The physical phenomena modeled include traveling and interacting shocks, acoustic waves, 
vortex shedding, viscous flows, detonation waves, cavitation, flows in fluid film bearings, heat conduction 
with melting and/or freezing, MHD vortex, hydraulic jump, crystal growth, chromatographic problems and 
solar wind [3-63]. In particular, its unexpected simple non-reflecting boundary conditions [12,50] and rather 
unique capability to resolve both strong shocks and small disturbances (e.g., acoustic waves) simultaneously 
[19,21,22] make the CESE method an effective tool for attacking computational aeroacoustics (CAA) prob- 
lems. Note that the fact that the second-order CESE method is capable of solving CAA problems accurately 
is an exception to the commonly-held belief that a second-order scheme is not adequate for modeling CAA 
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problems. Also note that, while numerical dissipation is needed for shock capturing, it may also result in 
annihilation of small disturbances. Thus a solver that can handle both strong shocks and small disturbances 
simultaneously must be able to overcome this difficulty. 

In spite of its nontraditional features and robust capabilities, the core ideas of the CESE method are 
simple. In fact, all of its key features are the results of a pursuit driven by these simple ideas. The first 
and foremost is the belief that the method must be solid in physics. As such, in the CESE development, 
conservation laws are enforced locally and globally in their natural space-time unity forms for ID, 2D and 3D 
cases. Moreover, because direct physical interaction generally occurs only among the immediate neighbors, 
use of the simplest stencil also becomes a common CESE feature. Obviously, this requirement has the effect 
of simplifying boundary-condition implementation. 

The second idea emerges from the realization that stability and accuracy are two competing issues in 
time-accurate computations, i.e., too much numerical dissipation will degrade accuracy while too little of 
it will cause instability. In other words, to meet both accuracy and stability requirements, computation 
must be performed away from the edge (“cliff”) of instability but not too far from it. This represents a 
real dilemma in numerical method development. As an example, high order schemes generally have higher 
accuracy and lower numerical dissipation. However, they are susceptible to computational instabilities. In 
fact, in complicated real-world applications, not only they seldom live up to their nominal order of accuracy — 
generally they possess only first-order accuracy when shocks are present, stability of these schemes often is 
difficult to maintain without resorting to problem-specific treatments. To confront this issue head-on, in 
the CESE development, generally it is required that a solver be built from a nondissipative (i.e., neutrally 
stable) core scheme. By definition, computations involving a neutrally stable scheme are performed right on 
the edge of instability and therefore the numerical results generated are nondissipative. As such numerical 
dissipation can be controlled effectively if the deviation of a solver from its nondissipative core scheme can 
be adjusted using some built-in parameters. 

Moreover, because an accurate viscous flow simulation requires that the numerical dissipation be much 
smaller than the physical dissipation which decreases as the Reynolds number increases, in principle, an 
accurate and robust solver for high-Reynolds-number Hows must be able to cut numerical dissipation as 
the Reynolds number increases. Obviously this requirement can only be met by a solver built from a 
nondissipative core solver. 

Other CESE ideas are: (i) the flux at an interface be evaluated in a simple and consistent manner; 

(ii) genuinely multidimensional schemes be built as simple, consistent and straightforward extensions of ID 
schemes so that multidimensional shocks can be captured without using any mesh alignment technique; 

(iii) triangular and tetrahedral meshes be used in 2D and 3D cases, respectively, so that the method is 
compatible with the simplest unstructured meshes and thus can be easily used to solve problems with 
complex geometries; and (iv) logical structures and approximation techniques used be as simple as possible, 
and special techniques that has only limited applicability and may cause undesirable side effects be avoided. 
Fortunately for the CESE development, realization of the above lesser ideas (i)-(iv) follows naturally from 
the first two core ideas. 

Building on its initial successes, efforts to refine and improve the CESE method have continued in the 
past few years. As an example, it was shown in [5] that the numerical dissipation of a dissipative extension 
of a CESE core scheme may increase to an intolerable level as the value of the CFL number decreases from 
its maximum stability bound. As such, in a case with a large CFL number disparity (e.g., a simulation 
with a highly non-uniform spatial mesh and a spatially independent time step), the performance sensitivity 
with respect to the CFL number can lead to a solution that is highly dissipative in a region where the 
local CFL number <C 1. Even though a remedy was suggested in [5], a simple and robust solution to this 
problem had not arrived until a family of new Courant number insensitive schemes [25,31,46-49,51,54,55,58] 
was developed with fresh insights. Note that these new schemes have one important advantage, i.e., all 
variables at each mesh point can be evaluated explicitly without resorting to solving a system 
of linear/nonlinear algebraic equations involving these local mesh variables even in applications 
where systems of multidimensional nonlinear PDEs are solved. 

As another example, even though they are accurate enough to solve CAA problems, second-order CESE 
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solvers are not capable of resolving fine flow structures within a boundary layer without using relatively fine 
meshes and/or meshes with large aspect ratios. To overcome this limitation, two neutrally stable 4th-order 
schemes, referred to as the a( 3) and o(4) schemes [59-61], were developed in 2006 and intended to serve as 
the core schemes of other high order CESE schemes. Unfortunately, the dissipative extensions of these new 
schemes turned out to have the same intractable stability problem that afflicts traditional high order schemes. 
As such it became obvious to the author that overcoming this problem requires a fresh new prospective. 

It will be shown in this paper that, with the aid of a conceptual leap, one can build from a given 2nd- 
order CESE scheme its 4th-, 6th-, 8th-,... order versions which have the same stability conditions of 
the 2nd-order scheme and also retain all other advantages of the latter scheme. Thus the high 
order CESE schemes so constructed would avoid the common shortcomings of traditional high order schemes 
including: (a) susceptibility to computational instabilities; (b) computational inefficiency due to their local 
implicit nature (i.e. , at each mesh point, need to solve a system of linear/nonlinear equations involving all 
the mesh variables associated with this mesh point); (c) use of large and elaborate stencils which complicates 
boundary treatments and also makes efficient parallel computing much harder; (d) difficulties in applications 
involving complex geometries; and (e) use of problem-specific techniques which are needed to overcome 
stability problems but often cause undesirable side effects. 

Because the new high order CESE schemes are built from 2nd-order CESE schemes, the remainder of 
the paper includes an extensive review of 2nd-order CESE schemes (Sec. 2-4). It is followed by a thorough 
discussion of several linear and nonlinear high order CESE schemes (Sec. 5-7). A sketch of multidimensional 
extensions is provided in Sec. 8. The numerical results are presented in Sec. 9. Conclusions and discussions 
are given in Sec. 10. As a final note, in Sec. 10, we also present a hybrid 3rd-order scheme built from the 
Lax-Wendroff scheme and a 2nd-order CESE scheme. 

2. Review of second-order ID CESE schemes 


For simplicity, we review existing 2nd-order CESE solvers for the simple partial differential equation 
(PDE) 

du du , „ , 

— + a— = 0 (2.1) 


dt 


dx 


where the advection speed a / 0 is a constant. Let X\ = x, and X 2 = t be the coordinates of a two- 
dimensional Euclidean space E 2 . Then, because Eq. (2.1) can be expressed as V • h = 0 with h = f ( au,u ), 
Gauss’ divergence theorem in the space-time E 2 implies that Eq. (2.1) is the differential form of the integral 
conservation law 

(f h-ds = 0 (2.2) 

JS(V ) 


As depicted in Fig. 1, here (i) S(V) is the boundary of an arbitrary space-time region V in E 2l and (ii) 
ds = dan with da and n, respectively, being the area and the unit outward normal vector of a surface element 
on S(V). Note that: (i) because h ■ ds is the space-time flux of h leaving the region V through the surface 
element ds, Eq. (2.2) simply states that the total space-time flux of h leaving V through S(V) vanishes; 
(ii) in E 2 , da is the length of a line segment on the simple closed curve 5(U); and (iii) all mathematical 
operations can be carried out as though E 2 were an ordinary two-dimensional Euclidean space. 

To proceed, let fi denote the set of all space-time staggered mesh points in E 2 (dots in Fig. 2(a)), 
where n = 0, ±1/2, ±1, ±3/2, ±2, . . ., and, for each n, j = n ± 1/2, n ± 3/2, n ± 5/2, . . .. Each (j, n) € O 
is associated with a solution element, i.e., SE(j, n). By definition, SE(j, n) is the interior of the space-time 
region bounded by a dashed curve depicted in Fig. 2(b). It includes a horizontal line segment, a vertical line 
segment, and their immediate neighborhood. By definition, the end points of the line segments referred to 
above are excluded from SE(j, n) so that two SEs will not overlap. 

Eq. (2.2) will be simulated numerically assuming that, for any (x,t) £ SE(j, n), u(x,t) and h(x,t), 
respectively, are approximated by 


* / , ■ \ def n 

u (x,t;j,n) = u 1 


± (Uxtfix - Xj) ± (u t )?(t - t n ) 


(2.3) 


3 



and 


h*(x,t]j,n) d =(au*(x,t-,j,n),u*(x,t-,j,n)) (x,t) £ SE (j,n) (2.4) 

Note that: (i) u", (u x )j, and (u t )" are constants in SE (j,n), and (in a rough sense) they can be considered 
as the numerical analogues of the values of u, du/dx, and du/dt at the mesh point ( j,n ), respectively, (ii) 
(. Xj,t n ) are the coordinates of the mesh point (j, n) with Xj = j Ax and t n = nAt, and (iii) Eq. (2.4) is the 
numerical analogue of the definition h = ( au, u ). 

Let u = u*(x,t;j,n ) satisfy Eq. (2.1) within SE Then one has 

(u t )? = -a(u x )? (2.5) 


As a result, Eq. (2.3) reduces to 


”{x,t; j, n) = u" + ( u x )j [(a; - Xj) - a(t - t n ) 


( 2 . 6 ) 


i.e. , Uj and (u x )" are the only independent mesh variables associated with (j, n ). Note that, for the purpose 
of uniquely defining the flux over a boundary of a CE, u*(x, t ; j, n ) is defined here only for (x, t) G SE (j, n ) 
for each ( j,n ) £ Cl. In some later developments, this function will also be used in applications other than 
flux evaluation such as construction of finite-difference approximations. Hereafter, for such extended 
applications, it is to be understood that u*(x,t',j,n .) is defined even for a point (x,t) </ SE (j,n). 

Let E 2 be divided into non-overlapping rectangular regions (see Fig. 2(a)) referred to as basic conser- 
vation elements (BCEs). As depicted in Figs. 2(c)-2(e), (i) each (j, n) € Cl is assigned with two BCEs, i.e., 
CE_(j, n) and CE + (j, n); (ii) each BCE has one and only one pair of diagonally opposite vertices which 
belong to O; (iii) the space-time E 2 can be filled by CE_(j, n) and CE + (j, n), ( j,n ) € f 2; and (iv) CE(j,n), 
which is the union of CE_(j, n) and CE + (j, n), is referred to as a compounded conservation element (CCE). 

Given the above preliminaries, we are ready to describe the existing CESE solvers for Eq. (2.1). 

2.1. The a scheme 

Note that, among the line segments forming the boundary of CE_(j, n), AB and AD belong to SE(j, n), 
while CB and CD belong to SE(j — 1/2, n — 1/2). Similarly, the boundary of CE + (j, n) belongs to either 
SE(j, n) or SE(j + 1/2, n — 1/2). As a result, by imposing the two local flux conservation conditions at each 
(j,n) G Cl, i.e., 

h*-ds = 0 (j,n)eCl (2.7) 

7s(CE + o») 


and 


h* ■ ds = 0 


s(CE_o») 


(. 3 . n) £ Cl 


(2.8) 


and using Eqs. (2.4) and (2.6), one can obtain two equations for two unknowns u" and (ita.-)^- fact, with 
the aid of 

v = a— and («*)? d A f (j, n) £ Cl (2.9) 

AX J 4 J 

it can be shown that (see Comment (a) given at the end of this subsection) Eqs. (2.7) and (2.8) are equivalent 
to 

1 ^ 2 (j,n)£Cl (2.7 a) 


(1 - v) [u + (1 + = (1 - v) [u - (1 + v)u x \ j+1 'i 2 


and 


- 1/2 


(1 + u)[u-( 1 - v)u x \ J = (1 + v) [u + (1 - 


(. j , n) G Cl 


(2.8 a) 


respectively. To simplify notation, in the above and hereafter we adopt a convention that can be explained 
using an expression on the left side of Eq. (2.7a) as an example, i.e., 


[it T (1 + v)u x L — u" + (1 + v){ux) T j 
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By adding Eqs. (2.7a) and (2.8a) together, one has 

= \ ( ^ ^ ^ + ^ u ^j+ 1/2 + (1 + [« + (1 - v)ux]™ll /2 } U, n) G O (2.10) 

Let 1 — v 2 yf 0, i.e., 1 — v yf 0 and 1 + v yf 0. Then Eqs. (2.7a) and (2.7b) can be divided by (1 — v) and 
(1 + v), respectively. By subtracting the resulting equations from each other, one has 

(«*)" = («*)" ( 2 . 11 ) 

with 

(«£)" = \ { [« - a + ^)^];; 1 1 / 2 2 -[«+ a - ^)^];: 1 1 / 2 2 } c ?, «) e n ( 2 . 12 ) 

It has been shown that Eqs. (2.7) and (2.8) imply (i) Eq. (2.10) for all and (ii) Eq. (2.11) if \v\ yf 1. 
On the other hand, it can be shown that Eqs. (2.10) and (2.11) are equivalent to 

[u + (1 + is)ux]j = [u - (1 + ^)mx]"+i/ 2 2 (i > n ) e ^ (2.76) 

and 

[u - (1 - v)ux]j = [u + (1 - (j, n) G O (2.86) 

for all v. Because (i) Eqs. (2.7b) and (2.8b) imply Eqs. (2.7a) and (2.8a), respectively; and (ii) Eqs. (2.7a) 
and (2.8a) are equivalent to Eqs. (2.7) and (2.8), respectively, one concludes that, for all v, the last two 
conservation conditions are enforced by the a scheme which, by definition, is formed by Eqs. (2.10) and 
(2.11) for all v. 

Note that, because 

du ax du 

dx 4 dx 

def 

if x = x/(ax/4), the normalized parameter (u s )" can be interpreted as the value at (j, n) of the derivative 
of u with respective to the normalized coordinate x. Also note that: (i) the superscript symbol “a” in the 
parameter (m|)” is introduced to remind the reader that Eq. (2.11) is valid for the a scheme, and (ii) the 
a scheme is second-order accurate in space and time for both u” and the non-normalized parameter (w x )” 

[i.5]- 

The review of the a scheme is concluded with the following remarks: 

(a) Evaluation of Eqs. (2.7) and (2.8) can be facilitated by the following observations: because u*(x, t ; j, n) is 
linear in x and t, the total flux of h* leaving CE_ (j, n) or CE + (J, n ) through any of the four line segments 
that form its boundary is equal to the scalar product of the vector h* evaluated at the midpoint of the 
line segment and the “surface” vector (i.e., the unit outward normal vector multiplied by the length) of 
the line segment. 

Note that, by definition, points B and D depicted in Fig. 2(c) do not belong to either SE (j,n) or 
SE(j — 1/2, n — 1/2) . This fact, however, does not pose a problem for flux evaluation over S (CE_ ( j , n)) 
because the values of h* at isolated points do not contribute to the flux of h* over a finite line segment. 
A similar comment applies to points D and F depicted in Fig. 2(d). 

(b) Even though they have completely different origins and structures (one is a two-level scheme with two 
independent variables per mesh point while another is a three-level scheme with a single independent 
variable per mesh point), the a scheme and the classical leapfrog scheme are deeply related. In fact, it 
can be shown that each of the two variables u)' and (u x )(j of the a scheme by itself satisfies the same 
leapfrog scheme automatically [56]. Moreover, after sorting through the differences in the notations 
used in typical CESE papers and traditional literatures, and also taking into account the decoupled 
nature of a leapfrog solution, it is shown in [5, pp.301 and 313-315] that the two amplification factors 
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of the a scheme are identical to those of the neutrally stable leapfrog scheme. Thus, the a scheme is 
also neutrally stable (i.e. , nondissipative) if \v\ < 1 (see the additional discussions given in Sec. 2.2). 

(c) Even though it is introduced to model a single PDE (i.e., Eq. (2.1)) with a single dependent variable 
u, the a scheme is formed by two coupled discrete equations (i.e., Eqs. (2.10) and (2.11)) involving 
two independent numerical variables u” and (u x )”. As such one would expect that the a scheme be 
consistent with a pair of PDEs with one of them being Eq. (2.1). Indeed, it is shown in [1] that, regardless 
how At > 0 and ax > 0 are refined, Eqs. (2.10) and (2.11) are always consistent with Eq. (2.1) and 
v — du/dx = 0 with the understanding that v being the analytical version of (u x )/j. 

(d) Because (i) the vector h* at any surface element lying on any interface separating two neighboring BCE 
is evaluated using the information from a single SE, and (ii) the unit outward normal vector on the 
surface element pointing outward from one of these two neighboring BCEs is exactly the negative of that 
pointing outward from another BCE, one concludes that the flux leaving one of these BCEs through the 
interface is the negative of that leaving another BCE through the same interface. As a result of this flux 
cancelation, the local flux conservation relations Eqs. (2.7) and (2.8) lead to a global flux conservation 
relation, i.e., the total flux of h* leaving the boundary of any space-time region that is the union of any 
combination of BCEs will also vanish. In particular, because CE (j,n) is the union of CE_(j, n) and 
CE +(j,n), 


s( CE(i») 


h* -ds = 0 


(j, n ) G O 


(2.13) 


must follow from Eqs. (2.7) and(2.8). In fact, it can be shown that Eq. (2.13) is equivalent to Eq. (2.10). 

(e) In addition to the nondissipative a scheme, as will be shown, there is a family of its dissipative extensions 
in which only the less stringent conservation condition Eq. (2.13) is assumed [5]. Because Eq. (2.13) is 
equivalent to Eq. (2.10), for each of these extensions, u” is still evaluated using Eq. (2.10) while (tii)" 
is evaluated using an equation different from Eq. (2.11). 

2.2. The a-e scheme and the c scheme 

Consider any ( j , n) G O. Then ( j ±1/2 , n — 1/2) G O. Let 


■j±i /2 = f u* (xj±y 2 ,t n ;j ± 1/2, n- 1/2) = [u + (At/2)u t \ 


- 1/2 
j± 1/2 


(2.14) 


The validity of the last equality sign in Eq. (2.14) follows from Eqs. (2.5) and (2.6). Thus u'j±i/ 2 is the 
first-order Taylor approximation at the point (xj^/ 2 ,t n ) (which ^ SE(j ± 1/2, n — 1/2)) evaluated using the 
mesh variables at ( j ±1/2, n — 1/2). Moreover, with the aid of Eqs. (2.6) and (2.9), one can recast Eq. (2.14) 
in the form 

Uj± 1/2 = {u-2v Wx)"±i/ 2 2 (2.15) 

Note that, by definition, (j ±1/2 , n) </ if (j,n) G f l. Thus w'^ 1 y 2 is associated with a mesh point ^ O. 
The reader is warned that similar situations may occur in the rest of this paper. 

Because v!?± 1 , 2 is a first-order Taylor’s approximation of u at ( j ± 1/2, n), 


(u%)7 = 


def AX I u j+ 1/2 U j- 1/2 


n ,n — ni' n 
“j+1/2 j — 1/2 


AX 


(2.16) 


is a central-difference approximation of du/dx at ( j,n ), normalized by the same factor ax/4 that appears 
in Eq. (2.9). Note that: (i) the superscript “c” is used to remind the reader of the central-difference origin 
of the term (u|)”; and (ii) by using Eqs. (2.15), (2.16) and (2.12), one has 


K)) 



2 vus) 


n-l/2 

J + 1/2 


(u 


2vu s ) 


n—l/2 
J — 1/2 


(2.17) 


6 



and 

K - <)] = | [(«+ 2 u^Zl',1 - (u - 2u,)]-$] (2.18) 

The a-e scheme is formed by Eq. (2.10) and 

(«*)" = («*)" + 2 c(«|-«|)7 (2.19) 

where e is an adjustable real number. Note that: (i) with the aid of Eqs. (2.12) and (2.18), Eq. (2.19) can 
be rewritten as 

(%)" = \ {[(! - e)« + (2e - 1 - - [(1 -e)u+( \-u- 2 e)u^r\^ } (2.19a) 

(ii) the a-e scheme reduces to the a scheme when e = 0; (iii) for the case e = 1/2, Eq. (2.19) reduces to 

(«*)? = («£)" ( 2 - 20 ) 

and (iv) because (w|)" represents a central-difference approximation, hereafter, to simplify its frequent 
references, the special a-e scheme with e = 1/2 will be referred to as the c scheme. 

As a preliminary to later developments, we offer the following remarks on the a-e scheme: 

(a) At each (j, n) G Q, Eqs. (2.10) and (2.11) imply Eqs. (2.7) and (2.8). Because Eq. (2.19) does not reduce 
to Eq. (2.11) except in the special case e = 0, at each (j, n) G f 2, generally the a-e scheme satisfies only 
the single conservation condition Eq. (2.13) (which is equivalent to Eq. (2.10)) rather than the two 
conservation conditions Eqs. (2.7) and (2.8). However, because (u|)™ generally is present on the right 
side of Eq. (2.19), implementation of the a-e scheme generally will still be burdened with the cost of 
solving two conservation conditions at each mesh point. The exception occurs only for the special case 
e = 1/2 (i.e., the c scheme) in which Eq. (2.19) reduces to Eq. (2.20). As it turns out, implementation of 
an Euler version of the c scheme does not require inverting any system of equations while that involving 
an Euler version of any other a-e scheme (e y^ 1/2) generally requires inverting, per mesh point and 
per time step, a system of several linear equations. As an example, because the 3D Euler equations is 
formed by a system of five PDEs with each of them being associated with three independent spatial 
derivatives at each mesh point, a 3D Euler a-e (e y^ 1/2) scheme requires inverting a system of 15 
(=5x3) equations. Partly because of this cost consideration, extensions of the c scheme have been 
used extensively. 

(b) For the a-e scheme, it is shown in [5] that the principal and spurious amplification factors per marching 
step (i.e., per At/2- recall that the solution is advanced by a time period At/2 per marching step) are 
A+ and A_, respectively, with 

A±(e, v, 9) = f ecos(0/2) — ivsm(0/2) ± (1 — e) [(1 — e)cos 2 (0/2) + (1 — i/ 2 )sin 2 (0/2)] (2.21) 

Here (i) i = f V^+T, and ( ii) 9, —oo < 8 < +oo, is the phase angle variation of a Fourier component over 
the length ax (i.e., 9 = kAx where k denotes the phase angle variation per unit length). In addition, it 
can be shown that the a-e scheme is stable <G> (i.e., “if and only if”) either 

0 < e < 1 and \v\ < 1 (2.22a) 

or 

e = 1 and \v\ = 1 (2.225) 

(Note: Because Eq. (2.11) is a result of Eqs. (2.7) and (2.8) only if \u\ y^ 1, the a and a-e schemes were 
not defined for the case \v\ = 1 in [5]. As such the case Eq. (2.22a) was not presented in [5]). Also 
one can show that the a-e scheme becomes more dissipative as the value of e increases from 0 to 1 [5]. 
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Note that, unless specified otherwise, in the remainder of the paper the ranges of e and v respectively, 
are defined by Eq. (2.22a,b). Moreover, note that: (i) the mesh points used in the a-e scheme are 
so arranged that those at any two time levels are staggered in space-time unless these two levels are 
separated by a whole integer multiple of At, and (ii) by applying Eqs. (2.10) and (2.19) successively 
over two consecutive marching steps, one can derive a compounded form (Eq. (3.7) in [5]) of the a-e 
scheme in which the mesh variables at any (j, n + 1) € are determined in terms of those at (j — 1, n), 
(j,n) and (j + 1 ,n) — i.e., the space-time staggering nature of the original scheme is eliminated from 
the compounded form. As such, (i) Hereafter two consecutive marching steps will be referred to as a 
marching cycle , and (ii) in the following accuracy evaluation of the numerical amplification factors, we 
choose to consider the amplification factors per marching cycle (i.e., per At). 

Note that, for several reasons which are related to the space-time staggered nature of the meshes used in 
the CESE method (see pp. 313-315 in [5]), the same spatial-interval length and time-step size denoted 
by Ax and At respectively in traditional literatures, are denoted by ax/2 and At/2 respectively here. 

• • d e f # a ' 

Because the same definition 6 = kAx is used here and in traditional literatures, for the same Fourier 
component (i.e., the same value of k), the corresponding value of 9 in traditional literatures is only one 
half of that in the current paper, i.e., the symbol 9/2 that appears on the right side of Eq. (2.21 ) would 
be replaced by 9 in traditional literatures. 

(c) Let k be a constant. Then u = e lk (. x ~ at ) represents a plane wave solution to Eq. (2.1). For this solution, 
k is the phase angle variation per unit length. Thus, by using the relations 9 = kAx and v = aAt/ ax, 
one concludes that 

e ik[x—a(t+At)] 

the analytical amplification factor per At = — jk(x-at ) — : e~ lkaAt = e~' lve (2.23) 

(d) The amplification factors of the a-e scheme per At (i.e., per two consecutive marching steps) are given 
by (A±) 2 . According to Eq. (2.21), [A±(0, v, 9)] 2 , the amplification factors of the a scheme (which 
corresponds to the case e = 0) per At, have the following properties: 



[A±(0, v, 9)} 2 \ = 1 if \v\<\ 

(2.24) 


hm [A+(0,zv, 9)] 2 = e Tl6 
12— *± 1 

(2.25) 


lim [A_(0, v, 9)] 2 = e ±i9 
12 — 1 

(2.26) 

and 

[A±(0, 0, 9)] 2 = 1 

(2.27) 

On the other hand, e“ 

lv9 , the analytical amplification factor per At, has the following properties: 


\e~ iv9 \ = 1 

(2.28) 


lim e~ iv9 = e Tiff 

v — 1 

(2.29) 

and 

e~ iv 6 = 1 if u = 0 

(2.30) 


For the a scheme, Eqs. (2.24)-(2.30) imply that: (i) the two amplification factor of the scheme, and 
the analytical amplification factor all have the same constant absolute value (= 1) if \u\ < 1, i.e., 
the scheme is nondissipative if \u\ < 1; (ii) in the limit of |^| — > 1 (i.e., v — > 1 or v — > —1), the 
principal amplification factor is identical to the analytical amplification factor and, thus, the former 
has no dissipative or dispersive error in this limit; (iii) also in the limit of \v\ — > 1, the phase angle 
associated with the spurious amplification factor is exactly the negative of that associated with the 
analytical amplification factor — thus the spurious amplification factor has a large dispersive error in 



this limit except when \6\ <C 1 (i.e., when the wavelengths of the errors 1); and (iv) when v = 0, the 
two amplification factors of the scheme, and the exact amplification factor are all equal to 1 and, thus, 
the two amplification factors of the scheme have no dissipative or dispersive error if v = 0. Because 
the accuracy of a scheme is essentially hinged on the behaviors of the principal amplification factor [1] , 
according to the facts stated above, the a scheme tends to become very accurate when \v\ approaches 
1 or 0. However, the dispersive errors associated with the nondissipative spurious amplification factor 
(which could be introduced at t = 0 as a result of an inaccurate initial- value specification [1]) may 
appear in a solution as persistent numerical wiggles when \v\ approaches 1. 

(e) According to Eq. (2.21), [A± (1/2, v, 9)] 2 , the amplification factors of the c scheme (which corresponds 
to the case e = 1/2) per At, have the following properties: 


2 _ x i0 


Y\m[\ + (l/2,v,e)\ = e 

v — 


lim [A_ (1/2, v, 6)\ 2 = — sin 2 (0/2) 


and 


[A±(l/2,0, 


2 L 


1 ± cos(0/2)\/2 — cos 2 (0/2) 


Because (i) e ln = —1 and (ii) 


1 - 


cos(0/2)\/2 — cos 2 (0/2) = [l — cos 2 (0/2)] >0 


(2.31) 

(2.32) 

(2.33) 


i.e., 


1 > 


cos(0/2)\/2 — cos 2 (0/2) 


Eqs. (2.32) and (2.33) can be recast as 


lim [A_ (1/2, v, 9)] 2 = sin 2 (0/2)e* 

I /— >±1 


(2.32a) 


and 


[A±(l/2, 0, 9)] 2 


1 

2 


1 ± cos (0/2) y/2 


cos 2 (0/2) 


(2.33a) 


respectively. By comparing Eqs. (2.29) and (2.30) with Eqs. (2.31), (2.32a) and (2.33a), one arrives at 
the following conclusions for the c scheme: (i) in the limit of \v\ — > 1, the principal amplification factor 
is identical to the analytical amplification factor and, thus, the former has no dissipative or dispersive 
errors in this limit; (ii) also in the limit of \v\ — » 1, the spurious amplification factor generally has large 
dissipative and dispersive errors; and (iii) when v = 0, both the principal and spurious amplification 
factors generally have large dissipative errors but no dispersive errors. According to the facts stated 
above, like the a scheme, the c scheme also tends to become very accurate when \v\ approaches 1. 
However, unlike the a scheme, the errors associated with the spurious amplification factor of the c 
scheme do die out rapidly when \ v\ approaches 1 if the value of\8\ is not too close to tt (see Eq. (2.32a)), 
i.e., if the wavelength of the discrete Fourier component is not close to 2ax — the shortest possible 
wavelength. Also, in sharp contrast to the a scheme, the c scheme becomes highly dissipative when v 
approaches 0. 

From the above discussions, one concludes that: 

(a) The advantages of the a scheme include: (i) it is nondissipative for all 9 if \v\ < 1; and (ii) when the value 
of \i/\ is close to 0 or 1, the scheme is very accurate. On the other hand, its disadvantages include: (i) 
because it is nondissipative, its extensions for nonlinear equations generally are unstable; (ii) when the 
value of 1 12 1 is close to 1, the dispersive errors associated with the nondissipative spurious amplification 
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factors will appear in a solution as persistent numerical wiggles; and (iii) comparing with the c scheme, 
it costs more to implement. 

(b) The advantages of the c scheme include: (i) when the value of \v\ is close to 1, it is very accurate and 
the errors associated with the spurious amplification factor generally die out rapidly; (ii) because of its 
dissipative property, its extensions for nonlinear equations can be stable; and (iii) in terms of ease of 
implementation and computer cost, it is much more superior than any other a-e scheme. On the other 
hand, the c scheme has a serious disadvantage, i.e., it is very dissipative when v approaches 0. 

Based on the remarks given above, it becomes clear that a new solver of Eq. (2.1) possesses all the 
advantages of the a and c schemes but none of their disadvantages if it is formed by Eq. (2.10) and a new 
equation in which (i) (u$)” is evaluated using a simple central-differencing procedure similar to that used 
to obtain (u|)"; and (ii) the (u s )" so obtained is identical to (w§)" in the limit of \v\ — * 1 and identical to 
( u % )” in the limit of v — > 0. In other words, such an ideal solver (i) is comparable to the c scheme in ease 
of implementation; (ii) becomes the c scheme in the limit of \v\ — > 1; and (iii) becomes the a scheme in the 
limit of v — > 0. In Sec. 3, it will be shown how such a solver can be constructed. 


2.3. The w-a scheme — a special wiggle-suppressing scheme 


If discontinuities are present in a numerical solution, any a-e scheme such as the c scheme is not equipped 
to suppress numerical wiggles that generally appear near these discontinuities. To serve as a preliminary 
for future development, here we shall briefly review an extension of the c scheme which was introduced as a 
remedy for this deficiency [5]. 

To proceed, let 


, def AX 

(«*-)" = — 


— li' n 

U j~ 1/2 


ii n — ti' n 

U 3 U J~l/2 


ax/2 


(2.34) 


and 

def AX ( u ' j '+ 1/2 ~ u j\ _ U 'j+ 1/2 “ 

[Us+) J ~ 4 ^ ax/2 J ~ 2 


(2.35) 


i.e., ( Ux _)" and (u®+)" are normalized numerical analogues of du/dx at (j, n) evaluated from the left and 
the right, respectively. By using Eqs. (2.16), (2.34), and (2.35), it can be shown that 


(«!)" = 2 («*- + «*+)" 


(2.36) 


i.e., (u§)" is the simple average of («$-)” and (u s+ )". As such, the c scheme can be extended by replacing 
(u|)" in Eq. (2.20) with an weighted average of (u$_)" and (u s+ )". In other words, the resulting extension 
is formed by Eq. (2.10) and 

(«*)? = («>-)"(«*-)" + («>+)?(«*+)? (2.37) 

where (■ w _)" and («;+)", the weight factors associated with ( u $ _)" and (u$+)” respectively, must satisfy the 
condition 

(«;_)" + («'+)" = ! ( 2 - 38 ) 

at all (j, n) £ O. In addition, it is assumed that the expression on the right side of Eq. (2.37) represents an 
interpolation (rather than an extrapolation) of (u s _)” and (u $+ )". This assumption <=> 

(w_)" > 0 and (w+) r - > 0 (2.39) 

For real variables x_, x+, and a > 0, let W- and W + be the functions defined by: (i) W_(x_, x + ; a) = 
W+(x_, x+; a) = 1/2 if x_ = x+ = 0; and (ii) 


W-(x-, x+; a) 


\x-\ a + \x+\ a 


(|x_| + |x+| > 0) 


(2.40) 
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and 


0-1 + 0 + | > o) 


(2.41) 


W+(x-, x+; a) 


|a’_| Q + 0+ 1“ 


if either ^ 0 or £+ ^ 0. Furthermore, let 


Oi)” = W ± ((ii*_)?, (u s+ )],a) (2.42) 

Then (u;_)" and (w+)" so defined satisfy Eqs. (2.38) and (2.39) and have the property that 

(«;_)" = («>+)" = 1/2 if a = 0 or | («*_)? | = |(u*+)?| (2.43) 

Note that: (i) to avoid dividing by zero, in practice either a small positive number such as 1CT 20 is added to 
each of the denominators in Eqs. (2.40) and (2.41), or one must make sure that the definition W-(x-, a) = 
W + (x-,x + ; a) = 1/2 instead of Eqs. (2.40) and (2.41) is applied when X- = x + = 0; and (ii) the special 
cases of Eqs. (2.40) and (2.41) with a = 1 and a = 2 are used in the slope-limiter proposed by van Leer [64] 
and van Albada et a 1. [65] . 

An extension of the c scheme is formed by Eqs. (2.10) and (2.37) with (w-)™ and (u>+)" being defined 
by Eq. (2.42). Because it involves an weighted average which is dependent on a parameter a, hereafter the 
scheme is referred to as the w-a scheme. In a straightforward manner, the w-a scheme has been extended 
to become the Euler w-a scheme. The details are given in [5]. 

Assuming a > 0 and |(u $ _)”| ^ |(« $+ )"|, then Eqs. (2.37), (2.40), and (2.41) imply that, of (w $ _)™ and 
(ux+)j‘, the one with smaller absolute value is assigned an weight factor > 1/2. This observation coupled 
with Eqs. (2.36)-(2.39) leads to the conclusion that, of («$-)” and («£+)", («$)” will have an algebraic 
value closer to the one with smaller absolute value if (u $ )™ is evaluated as an weighted average of (u £ _)” 
and ( u s +) 1 j according to Eq. (2.37). As a result, (u$)” so evaluated has a smaller absolute value than that 
evaluated using Eq. (2.20). In turn, numerical wiggles or overshoots can be annihilated by the additional 
numerical dissipation introduced as a result of this local “flattening” of (u$)"- It has been shown numerically 
that the extension is stable if \v\ < 1 and a > 0. Moreover, as a result of Eqs. (2.20), (2.36), (2.37) and 
(2.43), (i) the extension reduces to the c scheme when a = 0; and (ii) even if a > 0, the extension behaves 
very much like the c scheme in any smooth solution region (where the condition (u^-)!/ = (u$+)" more or 
less prevails) or at a solution extremum (where the condition (u$_)™ = ~(v,x+) , j more or less prevails). As 
such, the wiggle-suppressing power of the extension takes effect only if a > 0 and only in a solution region 
where |(u$_)!/| and | («£+)" | differ substantially. 

This section is concluded with the following discussion of CESE non-reflecting boundary conditions. 

2.4. Non-reflecting boundary conditions 

An unique feature of the CESE method is its simple and robust non-reflecting boundary conditions 
[12,50]. Let (j,n) be a mesh point on the right spatial boundary of the computational domain. Then (i) 
(j — 1/2, n— 1/2) represents an interior mesh point at the ( n — l/2)th time level; and (ii) 

u] = w"r 1/2 and (u $ )" = (2-44) 


form a set of non-reflecting boundary conditions. On the other hand, for a mesh point (j, n ) on the left 
boundary, (i) (j + 1/2, n — 1/2) represents an interior mesh point at the (n — l/2)th time level; and (ii) 


n n— 1/2 

U j ~ U j+ 1/2 


and (ux)j = (%)" + i/ 2 2 


(2.45) 


form a set of non-reflecting boundary conditions. The above non-reflecting conditions are applicable to all 
the second-order schemes described in Sec. 2-4. A rigorous justification of these and other simple CESE 
non-reflecting boundary conditions is given in [50] . 

3. The c-t and c-t* schemes 
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The ideal solvers of Eq. (2.1) defined at the end of Sec. 2.2 will be developed here. As a preliminary, we 
shall show that can be cast into a central-difference form when v = 0. 

According to Figs. 2(d), \EF\ — \AD\ = 0 when At = 0. Moreover, because of the assumption a/0, 
v = 0 -o- At = 0. Thus, in the case v = 0, Eq. (2.7) reduce to the condition that the sum of the fluxes leaving 
CE + (j, n) through its top face AF and bottom face ED is zero. Because (i) n = (0, 1) at AF and n = (0, —1) 
at ED ( n denotes the outward unit normal vector at any boundary of CE + (j, n)), (ii) \AF\ = \ED\, and 
(iii) h* = ( au*,u *), by using the flux evaluation rule referred to in Remark (a) given at the end of Sec. 2.1, 
one concludes that Eq. (2.7) reduces to 


(u + Ua 0 ” = (u- Ux)j + i /2 (v = 0 ) 

if v = 0. Similarly, Eq. (2.8) reduces to 

(u - Ux)™ = (u + (y = 0 ) 


(3.1) 


(3.2) 


if v = 0. Because Eqs. (2.10) and (2.11) (which form the a scheme) are equivalent to Eqs. (2.7) and (2.8) if 
\v\ Z 1 1, they must be equivalent to Eqs. (3.1) and (3.2) when v = 0. In fact, by subtracting Eq. (3.2) from 
Eq. (3.1), one obtains Eq. (2.11) with («/)” being the reduced form of Eq. (2.12) when v = 0, i.e., 


«) 


n 

3 



u-) n ~ 1/2 
x >j+ 1/2 


(u + u s ) 


n—1/2 
3 - 1/2 


( U = o) 


(3.3) 


On the other hand, by summing over Eqs. (3.1) and (3.2), one has the reduced form of Eq. (2.10) for the 
case v = 0. 

With the aid of Eq. (3.3) and the facts that: (i) (u — Ux)™^ represents the approximation of u at 

the midpoint of ED while (u + Ux)™Zy 2 represents that at the midpoint of CD (see Fig. 2(c,d)); and (ii) 
the distance between the two midpoints is ax/2, one concludes that, for the special case v = 0, (u§)" is 
indeed a central-difference approximation of du/dx at (j, n — 1/2) (which coincides with (j, n) when v = 0), 
normalized by the factor ax/4. QED. 

According to the above discussions, construction of the ideal solvers is hinged on finding central-difference 
approximations for (u 5 )” such that each approximation (i) becomes (w|)" in the limit of \u\ — » 1, and (ii) 
reduces to the expression on the right side of Eq. (3.3) when v = 0. As a result, these new solvers will be 
constructed as the subschemes of the c-r scheme, a new class of CESE solvers for Eq. (2.1) to be described 
immediately. 

3.1. The c-r scheme 

To proceed, refer to Fig. 3. Here M + and M~ denote the midpoints of AF and AB, respectively. 
Also P + and P~ are two points on BF that satisfy the following conditions: (i) P + is to the right (left) 
of A I + <£=> P~ is to the left (right) of M~\ and (ii) \M+P+\ = \M~P~~\. In addition, the parameter r is 
defined by: (i) r = 0 if P + coincides with M + \ (ii) r a.t/ 4 = |M+P+| > 0 if P + is to the right of M + \ and 

(iii) tax/4 = — \M+P+\ < 0 if P + is to the left of M + . Obviously, it follows from the above definitions 

that (i) r = 0 if P~ coincides with M - ; (ii) tax/ 4 = \M~P~\ > 0 if P~ is to the left of M~\ and (iii) 

t ax/4 = — \M~P~\ < 0 if P~ is to the right of M~ . 

Next let 

u(P+) d = u * ( x(P + ),t n ;j + 1/2, n- 1/2) (3.4) 

and 

u(P~) d = u* ( x{P~),t n ;j - l/2,n-l/2) (3.5) 

where x(P + ) and x(P~) denote the x-coordinates of P + and P ~ , respectively. Then, according to Eqs. (2.5) 
and (2.6), u(P + ) is the first-order Taylor’s approximation of u at P + evaluated using the mesh variables 
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at (J + 1/2, n — 1/2) (i.e., point E) while u(P ) is that at P evaluated using the mesh variables at 
(j — 1/2, n— 1/2) (i.e., point C). Moreover, with the aid of Eqs. (2.6) and (2.9) along with the relations 


X{P + ) = Xj + (1 + T)— = Xj+l/2 ~ (1 - t) — , 

and 

x(p+) = Xj - (1 + T )— = Xj_i /2 + (1 - r) — , 
Eqs. (3.4) and (3.5), respectively, can be simplified as 


— 00 < T < +00 


— 00 < T < +00 


u(P + ) = [u - (1 + 2is - T)ux] n j+ ljl 


(3.6) 

(3.7) 


(3.8) 


and 

u(P~) = [u + (1 - 2 V- t)us:}^ZI/2 (3.9) 

To proceed, note that: (i) the mesh point (j,n) (i.e., point A depicted in Fig. 3) is the midpoint of 
P-P+, and (ii) \P~P+\ = (1 + r) ax/2. Thus 


def AX 
(«*)" = — 


~ u(P + ) - u(P~) 
(1 + r ) ax/2 


_ u(P + ) - u(P~) 

2(1 + r) 




(3.10) 


represents a central-difference approximation of du/dx at the mesh point (j, n), normalized by the factor 
ax/4. Thus the new scheme defined by Eq. (2.10) and («$)" = (u s )” represents a solver for Eq. (2.1). 
Moreover, by using Eqs. (3.8)-(3.10), the second component of the new scheme can be explicitly expressed 
as 


(%)" = (Ux)j = 2(1 + T ) { l u - ( 1 + 2v - T ) U x\ n j+l /2 - [u+(l-2v- t)Ux]"_i /2 } ( T ± ^ 1 ) (3- 11 ) 

In view of the facts that (i) (u$)" represents a central-difference approximation of (%)", and (ii) the ap- 
proximation is associated with the parameter r, hereafter the new scheme will be referred to as the or 
scheme. 

Moreover, by using Eqs. (2.12), (2.18) and (3.11), one has 


(fix)? = (u%)) 


2 r 

1 + T 


(u% - um - 


v (i - t ) 

2(1 + r) L 


(Ux)j + i/2 - («x)?_i/2 I (t Y -1) 


n- 1/2 


(3.12) 


By comparing Eq. (3.12) with (2.19), one concludes that the c-t scheme generally is different from the a-e 
scheme. In fact, the or scheme becomes the a-e scheme <t=> either (i) r = 1 or (ii) v = 0. For the case r = 1, 
Eq. (3.12) implies that (u 5 )” = (u|)". In other words, the c scheme is the special case of the c-t scheme 
with r = 1, a fact that can also be deduced from the observation that the points P + and P~ depicted in 
Fig. 3, respectively, coincide with points F and B (i.e., the mesh points (j + 1/2, n) and ( j — 1/2, n)) if r = 1. 
On the other hand, when v = 0, the c-t scheme become the a-e scheme with e = 2r/(l + r). In fact one can 
further deduce that c-t scheme reduces to the a scheme <t=> v = r = 0. 

Because the c-t scheme is formed by two rather complicated equations (i.e., Eqs. (2.10) and (3.11)) 
involving two parameters v and r, it were not expected that its von Neumann stability conditions could be 
cast into an explicit analytical form. However, by using the Jordan canonical form theorem [66], it is shown 
rigorously in [46] that the c-t scheme is von Neumann stable 

v 1 < 1, t>t 0 (v 2 ), and (i/ 2 ,r)y^( 1,1) (3.13) 
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where 


if s = 0 


/ \ def 

To{s) = < 


4 - s - 2 v / 2(2 _ I - s 5 ) 


s-l + Vl-2s + 5s 2 
2s 


if 0 < s < rjj 


if — s — i 


Note that: 

(a) The current stability conditions agree completely with those generated through numerical experiments 
[31] (Note: the reader is warned that the condition r > T 0 (y 2 ) that appears in Eq. (3.13) here is 
expressed as r > r G ( \v\) in Eq. (3.12) of [31], i.e., the function r D defined here is different from the 
function r 0 introduced in [31]). 

(b) It can be shown that [46]: (i) r G (s) is continuous at s = 0, i.e., 

lim t 0 (s) = r 0 ( 0) = 0 

s— >0+ 

(ii) r 0 (s) is consistently defined at s = 3/11, i.e., both the second and third expressions on the right 
side of Eq. (3.14) implies that r 0 (3/ll) = 1/3; (iii) r 0 (s) is continuously differentiable at s = 3/11, i.e., 

lim r/(s) = lim r/(s) = 121/90 
s ^rr“ s ^rr + 

where t'(s) = f dr 0 (s)/ds; (iv) r 0 (s) is strictly monotonically increasing in the interval 0 < s < 1; (v) 
r 0 (l) = 1; and (vi) 

s < T 0 (s) < y/s, 0 < s < 1 (3.15) 

(c) For any given fixed value of |^| < 1, the c-t scheme tends to become more dissipative as the value of r 
increases from its minimum stability bound t 0 ( v 2 ). 

With the above preliminaries, the ideal solvers of Eq. (2.1) will be constructed in Sec. 3.2. 

3.2. The c - t * schemes 

The value of r used in the c-t scheme generally can be chosen independent of v. Here we will introduce 
a subset of the c-t scheme in which r is a function of v 2 for each member of this subset. As a preliminary, 
note that, by using the properties of r 0 (s) presented earlier, it can be shown that there exist infinitely many 
choices of a strictly monotonically increasing smooth function h(s ), 0 < s < 1, which possesses the following 
properties: (i) 

h{ 0) = 0; and lim h{s ) = 1 (3.16) 

S — > 1 — 


and (ii) 


h{s) > t 0 (s ) if 0 < s < 1 


By definition, the subset contains every special c-t scheme with 


r = h{v 2 ) 


(M<!) 


where h(s) is one of the special functions defined above. By using Eqs. (3.16)-(3.18) and the fact that h(s) 
is a strictly monotonically increasing function in the domain 0 < s < 0, one concludes that, for each member 
of the subset, (i) 

t = 0 if is = 0 (3.19) 


lim r = 1 
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(iii) r increases with \v\ in the domain 0 < \v\ < 1; and (iv) 

T > T 0 {v 2 ) (\v\ < 1) (3.21) 

By using the results presented in the above items (i)-(iii), one can easily infer from Fig. 3 a simple relation 
between the value of \v\ and the locations of P + and P~ , i.e., as the value of \v\ increases from 0 to 1, P + 
will move away from M + and edge toward the mesh point (j + 1/2, n) while P~ will move away from point 
M~ and edge toward the mesh point (j — 1/2, n). 

Recall that (i) (u s )" = (w|)" if r = v = 0; and (ii) (w$)" = (u|)" if r = 1. As such, Eqs. (3.19) and 
(3.20) imply that, for each member in the subset, (i) (u$)" = (u|)" if v = 0; and (ii) (u$)" = (u§)” in the 
limit of \i/\ — * 1 _ . In other words, in the domain \v\ < 1, all members in the subset are the ideal solvers 
defined at the end of Sec. 2.2. Moreover, by using Eq. (3.13) and (3.21), one can also show that these ideal 
solvers are also stable in the same domain. Hereafter, each of these ideal solvers will be referred to as a c-r* 
scheme. 

Corresponding to infinitely many choices of h, there are infinitely many different c-r* schemes. In 
particular, with the aid of Eq. (3.15), one can see easily that the special choice 

h(s) = y/s (0 < s < 1) (3.22) 

represents a strictly monotonically increasing smooth function which also satisfies Eqs. (3.16) and (3.17). 
Thus Eq. (3.18) implies that a c-r* scheme can be defined using the relation 

r=M (M<1) (3-23) 

Any c-r* scheme described above meets all the requirements of an ideal solver. However, because of 
stability problem, a c-r* scheme and its multidimensional extensions may not be robust enough for some 
complicated real-world applications. To overcome this difficulty, we introduce a special c-r scheme which is 
defined by Eqs. (2.10) and (3.11) with 

T = (3\v I (/3 > 1; M < 1) (3.24) 

where /? > 1 is a constant parameter. For this scheme, Eq. (3.19) is also valid. Thus it becomes the 
nondissipative a scheme when v = 0 and, therefore, it will not become overly dissipative as v — > 0 (i.e., the 
scheme’s numerical dissipation is CFL number insensitive). Moreover, because /? > 1, Eqs. (3.13) and (3.15) 
along with Comment (c) given following Eq. (3.15) imply that the scheme is stable if \v\ < 1 and it becomes 
more dissipative (i.e., more stable) as the value of (3 increases. However, also because (3 > 1, 

lim r = (3 > 1 (3.25) 

hhi- 


As such, except the special case (3 = 1, the scheme does not become the c scheme in the limit of \v\ — * 1 _ , 
i.e., it does not meet all the requirement of an ideal solver. 

Note that, unless specified otherwise, in the remainder of this paper we consider only the 
c-r* schemes or the special c-r scheme defined by Eq. (3.25). 

Next we will introduce weighted-averaging extensions of the c-r schemes along with two classes of 
advanced weighted-averaging techniques. 

4. Other c-r* extensions and new weighted-averaging techniques 


To proceed, let 

def ax ( v%-u(P-) \ = U™ - u(P~) 
x i 4 ^ (1 + t)ax/3J (1 + r) 


(4.1) 
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and 


(4.2) 


, def AX ( u(P+) ~ U™ \ _ U{P+) - u!? 

~~ 4 y (1 + t)ax/A J ~ (1 + r) 

Because |AP_| = \AP + \ = (1 + t)ax/ 4 (see Fig. 3), it is easy to see that (us-)™ and (m$+)” are two 
normalized one-sided difference approximations of du/dx at the mesh point (j,n) with one being evaluated 
from the left and another from the right. Also, it follows immediately from Eqs. (3.10), (4.1) and (4.2) that 

(fix)? = \ [(fix-)? + (&*+)?] (4-3) 

Moreover, by (i) substituting Eqs. (2.10), (3.8) and (3.9) into Eqs. (4.1) and (4.2), and (ii) using Eqs. (3.11) 
and (3.19), one arrives at the conclusion that 

(fix-)? = (fi*+)? = (fix)? iy = 0) (4.4) 


when v = 0. 

With the above preliminaries, several extensions of the c-r* scheme will be constructed in the following 
subsections. 

4.1. Scheme w-1 

A comparison of Eqs. (4.1)-(4.3) with Eqs. (2.34)-(2.36) reveals that an obvious extension of the c-r* 
scheme can be obtained by replacing (u s _)" and (u s+ )" in Eqs. (2.37) and (2.42) with (u s _)" and (« s+ )", 
respectively. In other words, the new extension is formed by Eq. (2.10) and 

(«*)? = (w-)?(fix-)? + («;+)?(«*+)? (4.5) 

with 

(w±)? = W±((fix-)?, (fix+)?, «) (4.6) 

Because the scheme is the first extension of the c-r* scheme in which (u $ )" is expressed as an weighted 
average of (iix-)™ and (ux+)”, for simplicity, hereafter it will be referred to as Scheme w-1. It has been 
shown numerically that Scheme w-1 is stable if \v\ < 1 and a > 0. However, as will be shown, Scheme w-1 
is deficient in one key aspect. 

Eqs. (4.4) and (4.6) along with the definition of W± imply that, for any given a > 0, (u>_)™ — > 1/2 
and (w+)” — > 1/2 as v — > 0. In other words, for Scheme w-1, the “weighted” average on the right side of 
Eq. (4.5) will approach the simple average as v — > 0. According to an explanation given in the last paragraph 
of Sec. 2.3, this implies that Scheme w-1 will lose its capability to suppress wiggles or overshoots when v 
becomes small. For this reason, even though the Euler version of Scheme w-1 performs much better than 
that of the w-a scheme in its ability to resolve shocks and contact discontinuities crisply in a wide range 
(from 1 to less than 0.001) of the global CFL number (i.e. , the maximal value of local CFL numbers), it is 
handicapped by the fact that wiggles or overshoots can appear near a discontinuity in a generated solution 
when the local CFL number there becomes extremely small. In the following, it will be shown that this 
weakness can be overcome by simple modifications of Eq. (4.6). 

4.2. Scheme w-2 

A new scheme, referred to as Scheme w-2 is formed by Eqs. (2.10) and (4.5) with (u>±)” being given by 
Eq. (2.42). In other words, although (u 5 )” is still constructed as an weighted average of (u s _)" and (% + )"j 
the associated weight factors (ty±)" are evaluated using (u s _)” and («$+)? (which are defined in Eqs. (2.34) 
and (2.35), respectively). Because the last two parameters, respectively, are identical to the special cases of 
(u s _)" and (■ Ux+) 1 j with r = 1 (see Eqs. (2.15), (2.34), (2.35), (3.8), (3.9), (4.1), and (4.2)), their values do 
not vary with v. As such, (w±)” 1/2 and therefore the weighted average on the right side of Eq. (4.5) will 

not turn into a simple average when v = 0. In other words, Scheme w-2 is still capable of annihilating the 
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numerical wiggles near a discontinuity even if v becomes small. It has been shown numerically that Scheme 
w-2 again is stable if \v\ < 1 and a > 0. 

Note that a possible drawback of Scheme w-2 is that the relation |(w $ _)"| < |(« $+ )”| (|(u 5 _)”| > 
|(u $+ )”|) does not automatically follow from |(u s _)"| < |(u $+ )”| (|(u s _)"| > |(u i+ )"|) and vice versa. 
As a result, at some local mesh points, it may happen that, of (%_)" and (%+)", the one with smaller 
absolute value may not be associated with a weight factor >1/2. According to a discussion given in the last 
paragraph of Sec. 2.3, this implies that there is no guarantee that, at all localities, the weighted-averaging 
induced numerical dissipation will be available to suppress wiggles or overshoots. Despite this possible failing, 
fortunately it has been demonstrated numerically that, not only are they capable of suppressing wiggles or 
overshoots robustly, Scheme w-2 and its Euler extensions are also highly accurate. 

In the following, schemes that overcome the weakness of Scheme w-1 and also avoid the theoretically 
possible failing associated with Scheme w-2 will be constructed using two new classes of weighted-averaging 
techniques more advanced than that given in Eqs. (2.40) and (2.41). 

4.3. New classes of weighted-averaging techniques 

To pave the way, first we shall discuss a limitation of Eqs. (2.40) and (2.41) as a generator of weight 
factors. 

Let x± ^ 0. Then, for a given a > 0, obviously W- —> 1/2 and W + — » 1/2 as \x + /x-\ — » 1. As such, 
when \x+/x-\ is very close to 1, then both W- and W + will be very close to 1/2 unless a >> 1. As a result, 
in case that (i) (us±)" 7^ 0, (ii) |(u s+ )"/(u S -)”| is very close to 1 (e.g., when v is very close to 0); and (iii) 
Eqs. (4.6) is assumed, then the only way to prevent the weighted average that appears on the right side 
of Eq. (4.5) from becoming almost a simple average (and thus causing the overshoot problem discussed in 
Sec. 4.1) is to increase the value of a used. However, this approach may be impracticable because numerical 
evaluation of a quantity such as x a for any real number x generally is hampered by round-off errors and 
thus becomes highly inaccurate if the value of a becomes too large, say 100. 

To overcome the limitation discussed above, two classes of weighted-averaging techniques much more 
potent and flexible than that discussed in Sec. 2.3 were presented in Sec. 4.3 of [31]. These new techniques 
can be used to constructed weighted averages of any set of N (> 1) parameters. By applying these new 
techniques for the special case N = 2 [62], one can construct two more extensions of the c-t* scheme. These 
extensions, referred to as Schemes w-3 and w-4, are presented in Sec. 4.4 and 4.5, respectively. 

4.4. Scheme w-3 


Let 


and 


«) 


def 



(Jo 

a= H 

O. 

o > 0 and v ^ 0) 

(4.7) 

!(«*+)? 

r - i(fix-)?i a 


(a > 0 and |(u s _)"| + |(ux+)"| > 0) 

(4.8) 

!(«*.+ )" 

i a + 1 (%-)" i i q 



n def 1 

< = 2 

[1 + 

(j{(j,n; «)] 

(4.9) 


f 

< 

if 0 < w™ < 1 



(*>'-)? = 

3 

1 

if w ? > 1 

(4.10) 


l 

0 

if w™ < 0 



«)" 

def ^ 

- («oy 

(4.11) 


where a 0 > 0 and a > 0 are preset parameters in the order of 1. Note that: (i) to avoid dividing by zero, in 
practice a small positive number such as 10~ 20 should be added to each of the denominators of the fractions 
that appear in Eqs. (4.7) and (4.8); and (ii) Eqs. (4.10) and (4.11) imply that 


(«/_)? + «)? = 1 


(4.12) 
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and 


(«/_)" — 0 and ( w +)7 — ® (4-13) 

i.e., (w'_)J and ( w represent a pair of interpolating weight factors. With the aid of these weight factors, 
Scheme w-3 is formed by Eq. (2.10) and 

(Ux)? = + «)"(«*+)" (4- 14 ) 

To explore Scheme w-3, note that an immediate result of Eqs. (4.7)-(4.9) is 

f > 1/2 if |(fix+)?| > | (&*-)?| 

= 1/2 if l(fi*+)"l = l(*3-)"l (^0) (4.15) 

l < 1/2 if |(«*+)7I < |(fi s _)?| 

In turn, by using Eq. (4.15), Eqs. (4.10) and (4.11) imply that 

f o < «)7 < 1/2 < k _)7 < i if |(u s+ )7i > |(«x_)7i 

4 (07 = (07 = 1/2 if |(&x+)7l = |(fix-)?| (* y^ 0) (4.16) 

( o < (07 < 1/2 < (07 < 1 if 1(0)71 < 1(0)71 

As such, of (ux ~) 7 and (ux+)™, the one with smaller absolute value is assigned an weight factor > 1/2 in 
Eq. (4.14). As explained at the end of Sec. 2.3, it follows that Scheme w-3 also has the capability to suppress 
wiggles or overshoots. In fact, as will be shown, it would retain this capability even when \v\ <C 1. 

According to Eqs. (4.4) and (4.8), (i) £(j,n;a) implicitly is a function of v, and (ii) £(j,n;a) — > 0 as 
v — > 0. Let ((j, n; a) varies linearly with \v\ in the immediate neighborhood of v = 0, i.e., 

lim (£0, n; a)/\v\) = £+ (j, n; a) and lim (£(j, n; a)/\v\) = CO, n; a) (4.17) 

v—>0+ v— > 0 “ 


where C (J, n\ a) and £ 0 (j, n: a) are parameters independent of v. 
Eq. (4.9) and using Eq. (4.17), one has 

lim w 1 } = -\l + (j,n\ot) 1 and lim w™ = 

i/— »o+ J 2 ' »o _ J 

By assumption, o 0 > 0. Thus, one concludes that 


Then, by substituting Eq. (4.7) into 
\ [l + ^Ui.n;a)] (4-18) 


lim 

M-o+ J 


l 


(4.19) 


if 


C O', n; a) 0 and (j, n; a) y^ 0 
In turn, with the aid of Eq. (4.19), Eqs. (4.10) and (4.11) imply that 


lim (w'_ , )7 y^ - and 
M-o+ v ^ 2 


lim (u/_)7^ 

M-o+ v 


(4.20) 


(4.21) 


if Eq. (4.20) is assumed. In other words, the expression on the right side of Eq. (4.14) generally would not 
reduce to the simple average of (v,x~) 7 and («x+)7> and as suc ^ > Scheme w-3 would retain its capability to 
suppress wiggles and overshoots even when \v\ <C 1. This conclusion was supported by the numerical results 
of Euler solvers [31] which were constructed using the weighted-average technique similar to that presented 
here. 

4.4. Scheme w-4 
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Scheme w-4 is originally presented in [25]. As will be shown in [62], the weight factors used in this 
scheme can be constructed using the second class of weighted-averaging techniques reported in Sec. 4.3 of 
[31]. Let 


(v±)7 = 


(Ux±) 


n la 
3 I 


min{|(u 5 ;._)"|“,|(u $+ )"| Q } 


- 1 


(a > 0 and \(ux-)j\ + |(w$+)?| > 0) 


and 


(*+)? = 

1 + ( 7 7-)? 

iy Y o) 

2 + a[(77 + )? + (??-)?] 

III 

1 

1 + f 0 7+)? 

iy Y o) 

2 + < r[(rj + )? + (r?_)?] 


(4.22) 

(4.23) 

(4.24) 


where a > 0 is a preset parameter in the order of 1 while a is that defined in Eq. (4.7). Note that: (i) to 
avoid dividing by zero, in practice a small positive number such as 10~ 20 should be added to the denominator 
of the fraction that appears in Eq. (4.22); and (ii) Eq. (4.22) implies that 


(v+Y 


C n-Yi 


> 0 if |(fix+)?| > |(fix-)?| 
= 0 if|(fix+)?|<|(fix-)?| 

(4.25) 

>0 if |(«x-)?| > |(fix+)?| 
= 0 if|(fix-)?|<|(fix+)?| 

(4.26) 

(*+)? + (fi>-)? = 1 

(4.27) 

1? > 0 and {w _)? > 0 

(4.28) 


and 


Eqs. (4.23)-(4.26) imply that: (i) 


and 


i.e. , («)+)" and (u)-)? represent a pair of interpolating weight factors; and (ii) 

0 < («>+)? < 1/2 < («>-)? < 1 if |(fix+)?| > I (fix-)? I 
(*+)? = (*-)? = 1/2 if |(fix+)?| = | (fix-)? I 

0 < («>_)? < 1/2 < («;+)? < 1 if |(us + )?| < | (fix-)? | 

Scheme w-4 is formed by Eq. (2.10) and 

(«x)? = (fii-)?(fix-)? + (fii + )?(fix + )? 




(4.29) 


(4.30) 


According to Eq. (4.29), of (u s _)? and (u s+ )” in Eq. (4.30), the one with the smaller absolute value is 
also assigned an weight factor > 1/2. As such, Scheme w-4 also has the capability to suppress wiggles or 
overshoots. In fact, if (r/±)” varies linearly with \v\ in the immediate neighborhood of v = 0, one would 
conclude from the argument that lead to Eq. (4.21) that in general 


lim >+)? y 2 \ 

|i/|— >0+ J 2 


and 


1 


lim (*-)? yf 

|i/|— >0+ J 2 


(4.31) 


As such the expression on the right side of Eq. (4.30) generally would not reduce to the simple average of 
(ux-)j and (u s+ )”, i.e., Scheme w-4 would also retain its capability to suppress wiggles and overshoots even 
when \v\ <C 1. This conclusion was supported by the numerical results of Euler solvers [31] which were 
constructed using the weighted-average technique similar to that presented here. 

A rather thorough review of several 2nd-order ID CESE schemes has been provided in Sec. 2-4. In the 
following, it will be shown that each of these second-order scheme can be “naturally” extended to become its 
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4th-order version. In Sec. 5, first we will construct the a-4 and a-e-4 schemes, i.e., the 4th-order extensions 
of the a and a-e schemes, respectively. 

5. The a-4 and a-e-4 schemes 


To construct the a-4 and a-e-4 schemes, again we consider Eq. (2.1). Also we will use (i) the same mesh 
structure and the same CEs and SEs defined in Sec. 2 and (ii) the same definitions of h*(x,t;j,n), v, and 
(u 5 )" defined in Eqs. (2.4) and (2.9). However, Eq. (2.3) is now replaced by 

u*(x,t;j,n ) d = Uj + (u x )j(x - Xj ) + ( - t n ) + ^{u xx )](x - Xj) 2 

+ (Uxt)j{x - Xj)(t - t n ) + ^{u t t)j(t - t n ) 2 + ^(u xx X )j(x - Xj ) 3 (5.1) 

+ 7 }(Uxxt)]{x - Xjf(t - t n ) + ^(u xtt )j(x - Xj)(t - t ”) 2 + ^{u ut )j{t - t n ) 3 

Here uf, (u x )(t, ( u t )j , (u xx )™, {u xt )j, (u tt )j {u xxx ) r j, ( u xxt )j , (u xtt )J, and (uut)J are constants in SE(j,n). 
They can be considered as the numerical analogues of the values of u, du/dx, du/dt, d 2 u/dx 2 , d 2 u/dxdt, 
d 2 u/dt 2 , d 3 u/d x 3 , d 3 u/dx 2 dt, d 3 u/ dxdt 2 , and d 3 u/dt 3 at the mesh point (j, n), respectively. As such 
u*(x,t; j,n) represents a 3rd-order Taylor’s approximation of u. 

For any (j,ri) G O, let u = u*(x,t;j,n ) be a solution to Eq. (2.1) in SE(j, n), i.e., 


du* (x, t ; j, n) i __ du* (x, t ; j, n ) _ n 
i a ~ — u 


(x,t) £ SE (j,ri) 


dt dx 

Then one has 

K)? = -a(u x )?, (u xt )] = -a{u xx )?, (■ u tt )] = a 2 (u xx )], 

('U’XXt) j ~ Q'iW'xxx) j 1 (V'Xtt) j = O' (V'XXX )?, {nut) j — ^ i^xxx) j 


(. j , n) G n 


(5.2) 


(5.3) 


Substituting Eq. (5.3) into Eq. (5.1), one obtains 


u*(x,t;j,n ) = u] + (u x )] [(x - Xj) -a(t- t n )] 

+ 7}( u xx)j [{x - Xj) — a (t — t n )] 2 + ^{u xxx )j [(x - xj) - a(t - t n )\ 3 


Thus, (i) u*(x, t ; j, n) is a 3rd-order polynomial function of [(x — Xj) — a(t — t n )] ; and (ii) u", (u x )J, (m X x)", 
and ( u xxx )" are the only independent mesh variables associated with (j,n). Note that, for the same reason 
presented immediately following Eq. (2.6), except for flux evaluation, Eq. (5.4) and Eq. (5.7) (see below) are 
applicable even for (x,t) ^ SE(j,n). 

Let 


* / . • \ uei 

u xx (x,t;j,n) = 


def d 2 u*(x, t ; j, n) 


dx 2 


(5.5) 


Then Eqs. (5.2) and (5.4), respectively, imply that 


du* x (x,t\j,n) du* xx (x, t; j, n) = 
dt Q dx 


(x,t) G SE (j,n) 


(5.6) 


and 

u *xx{ x ^h n ) = i u xx)] + (u xxx )j [(a;- xj) - a(t- t n )\ (x, t) G SE (j,n) (5.7) 

i.e., (i) u = u* xx (x,t\ j,n) is also a solution to Eq. (2.1) in SE(j, n); and (ii) u* x (x,t; j,n) is a first-order 
polynomial function of [(x — Xj) — a(t — t n )] . Note that the fact that Eq. (5.2) implies Eq. (5.6) is but the 
special case of the more general result that Eq. (2.1) implies 


d_ (d^u\ d (d 2 u\ 

dt \ dx 2 J a dx \ dx 2 J 
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(5.8) 



In fact, for a smooth function u(x,t), the conservation condition Eq. (2.2) implies the conservation condition 



• ds = 0 


(5.9) 


where 

r def d 2 h _ f d 2 u d 2 u\ 

xx dx 2 \ dx 2 ’ dx 2 J 


(5.10) 


As a result of the above observations on the “equivalence” relation between u and d 2 u/dx 2 , the reader 
should no longer be surprised by the fact that, by replacing the symbols u”, (w x )", u*(x,t; j,n), and 
h*(x,t;j,n) with (u X x)j'i (u xxx )^, u xx (x,t\ j,n), and h xx (x,t; j,n), respectively, Eqs. (2.4) and (2.6) would 
turn into 

Kx( x >t'J> n ) = ( a Kx( x ,td,n),u*,. x (x,t-,j,n)) (x,t) £ SE (j,n) (5.11) 


(i.e., the numerical analogue of Eq. (5.10)) and Eq. (5.7), respectively. This structural similarity strongly 
suggests that the mesh variables (u ra )" and (u XX xYj can be advanced in time by using one of the 2nd-order 
CESE schemes described earlier with the understanding that the role played byu™ , u*(x,t;j,n), and 

h*(x,t;j,n) in a 2nd-order CESE scheme be replaced by (uxxYj > i u xxx)^, u* x (x,t; j,ri), and h* x (x,t;j,n), 
respectively . Note that this simplified approach is justified by the observation that (u xx )(j and (u JXI )", 
respectively, appear in the 2nd-order and 3rd-order terms on the right side of Eq. (5.4) and therefore the 
4th-order accuracy requirement of u*(x,t; j,n) demands only that (u xx )j and (u xxx )(j, respectively, be 2nd- 
order and lst-order in accuracy, i.e., the same orders of accuracy achieved by u ” and (it x )" hi a 2nd-order 
CESE scheme. 

For the reason presented above, equations such as Eqs. (5.11) and (5.7) hereafter will be designated as 
the substitution forms of Eqs. (2.4) and (2.6), respectively. With the aid of this designation, in the following 
(uxx)j' and (uxxx)j' will be advanced in time using the a and a-e schemes. 

5.1. Time marching of (u xx )J and ( u xxx )j 

By imposing the substitution forms of Eqs. (2.7) and (2.8), i.e., 


® h* ■ ds = 0 

Js( CE +(J ») 

(j, n) eci 

(5.12) 

and 

C 



o 

II 

Tec 

"TS 

8 
* H 

T-C 

IT 

1 

H 

@»^ 

(j, n) G Cl 

(5.13) 

and using Eqs. (5.7) and (5.11), one can show that (i) 



{Uxx)‘j = 7} {(1- - v ) [ u xx - (1 + v)Uxxx]™+l /2 + (1 + u ) [ u xx 

+ (1 - ^)uxxx]"r i 1 / 2 2 } 

(j,n)Gfi (5.14) 

and, assuming \v\ ^ 1, (ii) 

(%«)" = (Kxx)j 

(j, n) G Cl 

(5.15) 

with 

in \ n def 4 f r / . \ iTl — 1/2 l 

\ u xxx)j = 2 _ (1 "1 v ) u xxx\jj r i /2 ~ 1 

Uxx (1 l / )'Uxxx]j_x/2 

(5.16) 


where v, (uj)", (it^)", and (uxxx)j are defined by Eq. (2.9) and 


(«**)? = f (^) 2 (u**)? and (u XX x)]= (^) 3 («***)? 0»Gfi (5.17) 
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It can be shown that Eqs. (5.14)-(5.15) are equivalent to the substitution form of Eqs. (2.10) and (2.11). 
As such, the time marching scheme formed by Eqs. (5.14) and (5.15) literally is the same a scheme defined 
earlier. Also, one can shown easily that Eq. (5.14) is equivalent to the substitution form of Eq. (2.13), i.e., 


h* ■ ds = 0 


s(CEo») 


(j, n) G O 


(5.18) 


Note that, as explained in Comment (d) given at the end of Sec. 2.1, Eq. (5.18) is a direct result of Eqs. (5.12) 
and (5.13). 

Next let 

ulx(x,t;j,n) d = (^) u* xx (x,t;j,n ) and (u xxt )(j, n) d = (■ u xxt )\ ■ (j, n) G Cl (5.19) 


V 4 J 


Then the substitution forms of Eqs. (2.14)-(2.18) can be expressed as 


and 


( u xx)j± 1/2 — u xx { x j±i/ 2 jt n \ j E 1/2, n 1/2) — [u xx + (At/2)u xx t\j ±1 j 2 

( U xx)j± 1/2 = ~ 2 V Uxxx) j ±1 j 2 

def AX f (^xx)j+ 1/2 (^xx)j — 1/2 \ (Axx)j-t- 1/2 i^xx) j — 1/2 


= — 

V 11 xxx ) j a 


AX 


( u c ) n — 2 (u 2vu — \ n ~ 1 / 2 _ ( u 2vu — 

\ a xxx)j — A l “XX * l/u XXX )j+i/ 2 \ a xx * t/ u xxx )j_ 1 i 2 


,n- 1/2 


4 1 


( u xxx u xxx)j — A ( u xx + 2u xxx )j_ 1 i 2 (u xx 2u xxx )j + 1 i 2 


\n- 1/2 


4 L 


(5.20) 

(5.21) 

(5.22) 

(5.23) 

(5.24) 


respectively. Note that the validity of the last equality sign in Eq. (5.20) can be proved using Eqs. (5.7), 
(5.17) and (5.19) along with the relation ( u xxt )" = —a(u xxx )'j which is part of Eq. (5.3). 

Moreover, let the symbol e be replaced by ei. Then the substitution form of Eq. (2.19) can be expressed 
as 

(Uxx*)? = Kxx)7 + 2ei(*4xx - <xx)j (5-25) 

Obviously Eq. (5.25) reduces to (i) ( u xxx )^ = (■ u% xx ) ) if ei = 0; and (ii) (rtxxx)j = (“lxx)j ^ £ i = 1/2- With 
the aid of Eqs. (5.16) and (5.24), Eq. (5.25) can be rewritten as 


(uxss)] = ^ {K 1 “ ei)«xx + (2ei - 1 - ^)uxxx]" + i/ 2 2 “ K 1 ~ ei)«xx + (1 - v - 2e 1 )uxxx]”_i/ 2 2 } (5.26) 

Obviously Eq. (5.26) is the substitution form of Eq. (2.19a) if the symbol e in the latter equation is replaced 
by ei- 

Because (i) (u xx ) r ) and ( u xxx )" can be explicitly evaluated in terms of the mesh variables at the (n — 
l/2)th time level using Eqs. (5.14) and (5.26), and (ii) the latter equations are equivalent to the substitution 
forms of Eqs. (2.10) and (2.19a), it has been shown that indeed (u xx )™ and (u xxx )™ can be advanced in time 
using the a-e scheme. 

This completes the description of the first step in the construction of the a-e-4 scheme. In the following, 
we describe the second step, i.e., the details of how it" and ( u x )" can be advanced in time. As a preliminary, 
we conclude this subsection with the following remarks: 

(a) As a result of its construction, the scheme formed by Eqs. (5.14) and (5.26) observes both the conser- 
vation relations Eqs. (5.12) and (5.13) if and only if e\ = 0. 
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(b) Because Eq. (5.14) is equivalent to Eq. (5.18), the scheme formed by Eq. (5.14) and (5.26) still observes 
the less stringent conservation relation Eq. (5.18) if ei y^ 0. 

(c) Because of how they are derived, by replacing the symbols ei, Uxx, and Uxxx with e, u, and respec- 
tively, Eqs. (5.14) and (5.26), respectively, will become Eq. (2.10) and (2.19a), the defining equations of 
the a-e scheme. 

5.2. Time marching of it” and (ux)™ 

Let 

(s + )j = u - (1 + v)ux H -Uxx o -Uxxx (j,n)GO (5.27) 

L d J ]■ 

(»-)” = f U + (1 - v)ux + — — U + U - Uxx + — iy )( 1 + v ) u ___ (j, n) G O (5.28) 

L 4 6 J ■ 

Then, assuming Eqs. (2.4), (2.9), (5.4), and (5.17), it can be shown that the conservation relations Eqs. (2.7) 
and (2.8) are equivalent to 

\fr , , \ 2(1 + v + v 2 ) (1 + v){l + v 2 ) ] ^n- 1 / 2 ] n n ; r, n \ 

(!-,)<[«+(! + v)Ux H ^ u S x H g Uxxx J - (s+) i+ i/ 2 j 0 (5.29) 

and 

(1 + Z')j U-( l-v)Ux + — g + V ^ Uxx - — + ^ ) Mxxs - (s-)"ri/ 2 2 | = 0 (5.30) 

respectively. 

By adding Eqs. (5.29) and (5.30) together, one has 

^ [(! - ^)(s+)”+i/2 + (! + "KOS/a] - !(«**)" & n) € o (5.31) 

Let 1 — zA y£ 0, i.e., 1 — ^ 0 and 1 + ^ yf 0. Then Eqs. (5.29) and (5.30) can be divided by (1 — v) and 

(1 + v), respectively. By subtracting the resulting equations from each other, one has 

(«*)? = K~ 4 )j (J,n)efi (5.32) 

where 

« 4 )? = \ [(«+);; i 1 / 2 2 - (O^Ti/a] - f («**)" - ^ (! + ^ 2 ) («***)? (J. «) G n (5-33) 

Hereafter the superscript “4” is used to distinguish the parameters introduced in the current construction 
of the 4th-orcler schemes from their counterparts that appear in the 2nd-order schemes. 

At this juncture, note that: 

(a) The zero-order and first-order terms in the third-order expansion on the right side of Eq. (5.4) are 
identical to those in the first-order expansion on the right side of Eq. (2.6), and (ii) Eqs. (5.31) and 
(5.32) were obtained by imposing the conservation relations Eqs. (2.7) and (2.8) and assuming Eqs. (2.4) 
and (5.4) while Eqs. (2.10) and (2.11) were obtained by imposing the same conservation relations and 
assuming Eqs. (2.4) and (2.6). As such, one would expect that Eqs. (5.31) and (5.32) differ from 
Eqs. (2.10) and (2.11) only in the appearance of the extra 2nd-order and 3rd-order terms in the former 
equations. With the aid of Eqs. (2.12), (5.27), (5.28), and (5.33), the reader can prove for himself that 
this is indeed true. 
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(b) Assuming Eqs. (2.4), (2.9), (5.4), (5.17), (5.27), and (5.28), it has been shown that Eqs. (2.7) and (2.8) 
imply (i) Eq. (5.31) for all is\ and (ii) Eq. (5.32) if \v\ ^ 1. On the other hand, it can be shown that 
Eqs. (5.31) and (5.32) are equivalent to 


, 2(1 + is + is 2 ) (: l + v )(l + v 2 ) 

W 4 (1 d is)u x 0 Uxx 0 l^xxx 


/ \n— 1/2 

( s +)j + l/2 — 


(5.29a) 


J 3 


and 


■ - (1 - v)Ux + 


2(1 — is + is 2 ) 


(1-I/X1+!/ 2 ) 


- = 0 


J 0 


(5.30a) 


for all v. Because (i) Eqs. (5.29a) and (5.30a) imply Eqs. (5.29) and (5.30), respectively; and (ii) 
Eqs. (5.29) and (5.30) are equivalent to Eqs. (2.7) and (2.8), respectively, one concludes that, for all is, 
the last two conservation conditions are enforced by Eqs. (5.31) and (5.32). 

(c) Under the current assumptions, it can be shown that Eq. (5.31) is equivalent to the conservation relation 
Eq. (2.13). 

For Eqs. (5.31) and (5.32) to be considered as a complete marching scheme for «” and (%)", the 
mesh variables (wxx)j and (uxxx)j which appear on right sides of Eqs. (5.31) and (5.33) must be expressed 
explicitly in terms of the mesh variables at the (n — l/2)th time level. In the current construction, this is 
done by demanding that the former mesh variables be evaluated in terms of the latter mesh variables using 
Eqs. (5.14) and (5.26). By using this assumption and Eqs. (5.27) and (5.28), Eqs. (5.31)-(5.33) can now be 
cast into the forms by which and (it $ )" can be evaluated in terms of the mesh variables at the (n — l/2)th 
time level, i.e., 


l-is 


, \ 2 i/(l + is) 

U~( 1 + ls)Ux H 3 Usx 


(l + 1 s)(l-u 2 )^ 


1 - 1/2 


-I 3 + 1/2 


n- 1/2 
3 ~ 1/2 

and 

(u s )]=[u°r 4 (e i)]; (J»efl 

where 


l + is 


„ , 2is(l-is) (l~v)(l-v z ) 

U + (1 is)u x Uxx t , Uxxx 


(j, n) G n 


(5.34) 


(5.35) 


[«r 4 ( e i )]” d = - (! + v ) u x + 


1 + 3z/ 2 + ei(l + is 2 ) 


2\v(l — is 2 ) — ei(l + is 2 ) 

0 V'XXX 


n— 1/2 

J j+1/2 


1 

2 


U + (1 — is)ux + 


1 + 3i/ 2 + ei(l + is 2 ) 


2\u(l — is 2 ) + ei(l + is 2 ) 

l^XXX 


>. n- 1/2 
J 3 - 1/2 


(j, n) en 


(5.36) 

At this juncture, note that: 

(a) As a result of how Eqs. (5.34)-(5.36) are derived, one concludes that: (i) Eqs. (5.14) and (5.34) satisfy 
Eq. (5.31) automatically; and (ii) Eqs. (5.14), (5.26), and (5.35) satisfy Eqs. (5.32) automatically. 
Because the conservation relations Eq. (2.7) and (2.8), in turn, are enforced by Eqs. (5.31) and (5.32) 
for all is (see Comment (b) given following Eq. (5.33)), the two conservation relations are observed by 
the scheme formed by Eqs. (5.14), (5.26), (5.34), and (5.35) for all is and ei. 

(b) It is pointed out earlier (see comments (a) and (b) given at the end of Sec. 5.1) that the scheme formed 
by Eqs. (5.14) and (5.26) observes (i) the conservation relations Eqs. (5.12) and (5.13) if = 0; and 
(ii) only the conservation relation Eq. (5.18) if ei 0. 


24 



(c) As a result of the comments made in the above items (a) and (b), one concludes that the special scheme 
formed by Eqs. (5.14), (5.26), (5.34), and (5.35) with ei = 0 observes the four conservation relations 
Eqs. (2.7), (2.8), (5.12), and (5.13). In fact, for this special scheme, the time evolution of the four 
independent mesh variables u ”, (%)", («$$)", and ( u^xx )” is completely determined by these four 
conservation relations. As such, this special scheme is referred to as the a-4 scheme. 

(d) On the other hand, for the case e\ / 0, the scheme defined in the above item (a) satisfies only three 
conservation relations, i.e., Eqs. (2.7), (2.8), and (5.18). In the following, this scheme will be extended 
so that the extension satisfies only the two least stringent conservation relations Eq. (2.13) and (2.18). 

To proceed, first we will construct the current version of («§)" which was defined in Eq. (2.16). As a 
preliminary, Taylor’s expansion is invoked to obtain 


u(xj ± [(1 + t)ax/4], t n ) = jl ± [(1 + t)ax/4]c>/<9x + (1/2) [(1 + t)ax/ 4\ 2 8 2 /dx 2 

(r > -1) (5.37) 

± (1/6) [(1 + r) Ax/4] 3 d 3 /dx 3 + (1/24)[(1 + t)ax / 4} A 8 A / dx A } u( Xj , t n ) + 0[(ax) 5 ] 

Here (i) r > — 1 is an adjustable parameter, and (ii) (. Xj + [(1 + t)ax/4], t n ) and (xj — [(1 + t)ax/4], t n ), are 
the coordinates of points P + and P~ depicted in Fig. 3, respectively. By subtracting the two equations in 
Eq. (5.37) from each other and then rearranging the resulting equation, one has 

du ^ _ u(xj + [(l+T)Ax/4] 1 t n ) -u(x 3 - [(1 + t)ax/ 4], t n ) [(1 + t)ax] 2 d 3 u ^ r, , 4 i 

dx [ } (1 + t)ax/2 96 dx^ 3 ' )+ ^ 

(5.38) 

Next, by using Eqs. (3.6) and (3.7), one has 

u(xj ± [(1 + T)Ax/4],t n ) = u(xj- tjL/2 T [(1 — t)ax/4], + At/2) (5.39) 

The expression on the right side of Eq. (5.39) can be expressed as a Taylor’s expansion about the point 
( x j±i/ 2 ,t n ~ 1 ' 2 ). Substituting the expansion into Eq. (5.39), one has 

u ( Xj ± [(1 + t)ax/4], t n ) = u i ( P ± ) + 4th-order terms involving At and Ax (5.40) 


where 


u 4 (p±) d| f |x _|_ |(At/2)9/5t =F [(1 — t)ax/ 4]<9/(92’} + (l/2){(At/2 )d/dt =F [(1 — r)Ax/4\d/dxY 
+ (l/6){(Af/2)<9/<9f =F [(1 - r)Ax/4]d/dx} 3 |u(x j±1/2 ,t" _1/2 ) 


(5.41) 


Hereafter, we consider only a mesh refinement procedure in which the ratio At/ Ax and the parameter 
t are held constant as At — > 0 and ax —> 0. Because the advection speed a is a constant, in such a mesh 
refinement procedure, v is also held constant as At — > 0 and ax — > 0. Moreover, with this assumption, the 
4th-order terms which appear on the right side of Eq. (5.40) can be expressed as O [(ax) 4 ] , i.e., Eq. (5.40) 
can be expressed as 

u ( Xj ± [(1 + r)Ax/4],t n ) = u 4 (P ± ) + O [(ax) 4 ] (5-42) 

In turn, by substituting Eqs. (5.42) into Eq. (5.38), one has 


du _ u 4 (P+) ~u 4 (P ) 

dx 3 ’ (1 + t)ax/2 


[(1 + t)ax] 2 d 3 u 


96 


dx 3 


(xj,t n ) + O [(ax) 3 


(5.43) 
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For any (j,n) £ Cl, u 4 (P + ) and u 4 (P ) will be approximated by 


u 4 *(P+) = u*{ Xj + [(1 + r)Ax/4],f"; j + 1/2, n- 1/2) 

= [u - (1 - T + 2u)ux + (1/2) (1 - r + 2^) 2 - (1/6)(1 - r + 2A) 3 ” /+ ^ 

u 4 *{P~) d =u*{xj - [(1 + t)ax/ 4], f"; / — 1/2, n — 1/2) 

= [u + (1 - r - 2 v)u s + (1/2) (1 - r - 2v) 2 u S x + (1/ 6) (1 - r - 2A) 3 


respectively. According to Eqs. (3.6), (3.7), and (5.1)-(5.4), u 4 *(P + ) is the 3rd-order Taylor’s approximation 
(i.e. , the error is O [(ax) 4 ] ) of u at P + evaluated using the mesh variables at (j + 1/2, n — 1/2) while u 4 * ( P ~ ) 
is that at P~ evaluated using the mesh variables at (j — 1/2, n — 1/2). Note that the validity of the last 
equality signs on the right sides of Eqs. (5.44) and (5.45) follows from Eqs. (2.9), (3.6), (3.7), (5.4) and 
(5.17). Moreover, according to Eq. (5.43), the parameter 


n def U 4 *(P+)-U 4 *{P~) 
l (1 + t)ax/2 


[(1 + r) ax] 2 
96 


(u X xx 


h 


(5.46) 


represents a numerical analogue of the value of du/dx at (j, n) with 3rd-order accuracy. This is consistent 
with the 4th-order accuracy requirement of u*{x,t\ j,n) because (u x )/j (the numerical analogue of du/dx) 
appears in the first-order term on the right side of Eq. (5.4). By using Eq. (5.17) and 


T 4 / \l n def AX r a / \ 

[Uxivij = — [u x {r) 


(5.47) 


Eq. (5.46) can be cast into its normalized form, i.e., 


n def U 4 *(P+) - U 4 *{P~) 

J ~~ 2(1 + r) 


(1 + r) 




' ('Uxxx) i 


(5.48) 


Because points P + and P coincide with the mesh points (j + 1/2, n) and (j — 1/2, n), respectively, 
when r = 1, the current version of ( u |)™ is defined by 


/ c-4\n def 
\ a x )j ~ 


< 4 ( 4 )] 


n 

j 


1 r „ o 4 ^ 3 

-[u~2uux + 2u Uxx-—Uxxx 


-| n- 1/2 
- 3 + 1/2 


1 

4 L 


u — 2vux + 2v z Uxx ~ 


4 V 3 ] n ~ 1/2 _l( x „ 

^ xxx Yj*} Q\V J xxx)j 


(j,n)£Cl (5.49) 


The validity of the last equality sign in Eq. (5.49) can be proved by using Eqs. (5.44), (5.45), and (5.48) with 
t = 1. Substituting Eq. (5.26) into the last term on the right side of Eq. (5.49) and denoting the resulting 
expression as [it§" 4 (ei)] , one concludes that 


w $’ 4 ( £ i )]” = j\ u -2vux + 


— — u — 2v Ux + 


o .,2 4 ( 1 - e i)L. , 4 (1 + v — v 3 — 2ei) _ 

'U'XX I r» ^XXX 


2v z - 


4(1 — e x )i 4(1 - v + v* - 2ei) 


Uxx 


n-1/2 

3 + 1/2 
n-1/2 

3 - 1/2 


(j,n) £ Cl (5.50) 
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By definition, the a-e-4 scheme is formed by Eqs. (5.14), (5.26), (5.34), and 

(%)" = [«s' 4 ( e i)] j + 2e 2 [wj 4 (e l) - u a £ 4 {e i)]” (j» G ft (5.51) 

where e 2 is another adjustable parameter. With the aid of Eqs. (5.36) and (5.50), Eq. (5.51) can be expressed 
as 


( u x)J — — ^3(1 — 62)11 — 3(1 + ^ — 2e2 )tt$ + [l + 3v 2 — 6e 2 + (1 + ^ 2 )ei + 2(1 — ^ 2 )eiC 2 ] Uxx 

N n- 1/2 

+ 2[z/(l — u 2 ) + 2 e 2 — (1 + ^ 2 )ei — 2(1 — v 2 )e\e 2 \uxxx r 

J j+l/2 

(5.52) 

— — ^ 3(1 — C 2 )u + 3(1 — v — 2 e 2 )us + [l + 3 v 2 — 6 e 2 + (1 + ^ 2 )ei + 2(1 — i' 2 )ci£ 2 \uxx 

n n- 1/2 

+ 2[i/(l — v 2 ) — 2e2 + (1 + ^ 2 )ei + 2(1 — v 2 )ei£ 2 \uxxx r (j> n) G ft 

J 3 -l/2 


We conclude this subsection with the following remarks: 

(a) Eqs. (5.51) reduces to (i) (u s )" = [u%~ 4 (e i)]" if C 2 = 0; and (ii) (u$)" = [u§' 4 (ei)]" if e 2 = 1/2. Because 
the a- 4 scheme is the special scheme formed by Eqs. (5.14), (5.26), (5.34) and (5.35) with ei = 0, the 
a-e-4 scheme reduces to the a - 4 scheme when e\ = 62 = 0. 

(b) Because of how they are derived, by replacing the symbol e 2 with e and dropping all the 2nd-order and 
3rd-order terms, Eqs. (5.34) and (5.52), respectively, will become Eqs. (2.10) and (2.19a), the defining 
equations of the a-e scheme. 

In the following, we will discuss the stability of the a-e-4 scheme. As a preliminary, the scheme will be 
cast into a matrix form. 

5.3. Matrix form of the a-e-4 scheme 

The component equations of the a-e-4 scheme, i.e., Eqs. (5.14), (5.26), (5.34) and (5.52), can be cast 
into the matrix equation 


q{j,n) = Q-(ei, e 2 , v) q(j + 1/2, n- 1/2) + Q+(e 1 ,e 2 ,i') q(j - 1/2, n - 1/2) (j,n) G ft (5.53) 


where 


and 


with 


/ \ 

(«*)? 

(Uxx)j 


V {Uxxx)j / 


(j, n) G ft 


Q±i e li e 2 


/ qti 

4 

i 

i 

921 

922 

0 

0 

V 0 

0 


4 

4\ 

± 

i 

923 

924 

i 

i 

933 

934 

4 

9^4 ' 


+ def 1 ± V + def . 1 - V 2 i def v( 1 - V 2 ) ± def (1 - V 2 ) 2 

qfi = ~2~ > 9i 2 = ±—3—. qf : 3 = t — 3 — , qfi = t — g — 


(5.54) 


(5.55) 


(5.56) 
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± def 1 — £2 ± def 2e 2 — 1 A V _j_ def 1 + 3zA — 6e 2 + (1 + Z^ 2 )e 1 + 2(1 — Z' 2 )ei£ 2 

921 — “F ^ > 922 — ^ i 923 = "F g 

± def Kl - ^ 2 ) -F 2e 2 ± (1 + z/ 2 )ei ± 2(1 - ^ 2 )eie 2 

9 24 - "F 3 

i def 1 ± V ± def ,1 ^ V 2 ± def 1 ~ , def 2ei — 1 ± V 

933 “ 2 ’ 934 — ± 2 ’ 943 _ 2 ’ 944 _ 2 


(5.57) 


(5.58) 


For simplicity, hereafter (2 + (ei, e 2 , z') and Q_(ei, e 2 , z') may be abbreviated as Q+ and Q_, respectively. 

Eq. (5.53) represents the form of the a-e-4 scheme by which mesh variables can be advanced over At/2. 
By applying Eq. (5.53) successively over two consecutive marching steps, one obtain a compounded form by 
which mesh variables can be advanced over At, i.e., 


q(j, n+ 1) = {Q+) 2 q(j - 1, n) + (Q+Q- + Q-Q+)q(j, n) + ( Q-) 2 q{j + 1, n) {j, n) G O (5.59) 


Eq. (5.59) has the advantage that the mesh variables involved are associated with mesh points which are not 
staggered in space-time. 

To further explore the relations between the a-e-4 and a-e schemes, note that the component equations 
of the a-e scheme, i.e., Eqs. (2.10) and (2.19a), can be cast into the matrix form 

q(j, n ) = M_(e, v) q(j + 1/2, n - 1/2) + M+(e, v) q(j - 1/2, n - 1/2) (j, n) G O (5.60) 


where 

g(j, n) d = ( ] (j, n) G Q (5.61) 

V («*)?/ 

and 

/ l±v ±(l-i/ 2 )\ 

M±(e, v) = f - ( ] — oo < e< +00 (5.62) 

\ T(1 — e) 2e-l±z ,) 

With the aid of Eqs. (5.56)-(5.58) and (5.62), it can be shown that 

= M±(e 2 , v) (5.63) 


= M±(e i.i/) (5.64) 

By comparing Eqs. (5.63) and (5.64) with Eq. (5.55), one now arrives at the conclusion that the 2x2 
submatrices in the top-left and bottom right corners of the 4x4 matrix Q±(ei, e 2 , z') are M±(e 2 ,^) and 
M±(e 1 , zx), respectively. 

Because of the relations Eqs. (5.63) and (5.64), as will be shown in the next subsection, the a-e-4 and 
a-e schemes share essentially the same stability conditions. 

5.4. Stability conditions of the a-e-4 scheme 

Note that a scheme such as Eq. (5.53) would become computationally unstable if the round-off errors 
introduced initially at any time level are amplified without bound as they propagate down the subsequent 
time levels. To be specific, consider the effect of contamination of the initial values by round-off errors. Let 
q(j, n ) denote the solution associated with the given exact initial values while q(j, n) + 6q(j, n) denote the 
solution associated with the contaminated initial values. Then q(j,n .) must satisfy Eq. (5.53). Also, because 


and 


9u 9j t 2 


, 921 922 , 


9*3 9*4 


,943 ot. 
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the coefficient matrices Q+ and Q_ in the linear scheme Eq. (5.53) have constant elements independent of 
mesh variables, the contaminated solution satisfies the equation 


q(j, n) + 6q(j, n) = Q _ [q{j + 1 / 2 , n - 1 / 2 ) + Sq{j + 1 / 2 , n - 1 / 2 )] 
+ Q+ [q{j - 1 / 2 , n - 1 / 2 ) + Sq(j - 1 / 2 , n - 1 / 2 )] 

By subtracting Eq. (5.53) from (5.65), one has 

Sq(j, n) = Sq(j + 1/2, n - 1/2) + Q + 6q(j - 1/2, n - 1/2) 


(j, n) G O 


(j, n ) G ft 


(5.65) 


(5.66) 


Thus the time evolution of Sq(j, n), which represents the solution errors induced by the initial- value contam- 
ination, is also governed by the scheme Eq. (5.53). For this reason, stability of the linear scheme Eq. (5.53) 
can be determined by the time-evolution behavior of q{j, n) itself. 

The stability of the a-e-4 scheme will be studied using the von Neumann method. In this method (see 
Sec. 4 in [1] for a rigorous treatment based on the discrete Fourier analysis), the time-evolution behavior of 
each Fourier component (characterized by a phase angle 0) of q(j, n) is studied separately. Because these 
components propagate in time independently for a linear scheme such as Eq. (5.53), the scheme is stable if 
only if each component is bounded as it evolves in time. 

Let 


Q(J,n) = q*(n, 9)e 


ijO 


( j , n) G fl; —oo < 6 < +00; i = f \f—l 


(5.67) 


Here (i) 9 is the phase angle variation of a Fourier component over the length Ax, and (ii) q* (n, 9) is a 2 x 1 
column matrix. Substituting Eq. (5.67) into Eq. (5.59), one has 


q*{n + 1 , 9) = [Q(ei, e 2 , v, 6)] q*(n,9) 


n = 0,±l/2,±l, ±3/2, . . . 


(5.68) 


where 


Q(e 1 ,e 2 ,v,9) 

= e i6 / 2 Q-(e\, 

£ 2 ,^) + e i 0 // 2 Q+(e 1 . 

,Z 2 ,v) 




/e* el2 q : n 

+ e~ ie / 2 q+ < 

jei2 q - 2 + e -iei2 q + 2 

e iOI2q- 3 + e -iei2 qt3 

e iei 2 q lA + e -ie/ 2 q+ 4 

\ 


e i9 ' 2 q 2 1 

+ e~ ie / 2 q+ , 

js/2 q - 2 + e -^ 2 q + 2 

ei S / 2q- 3 + e -^ / 2q+ 

e i6 l 2 q24 + e_i6l / 2( ?24 




0 

0 

e ie/2q- 3 + e -ie/2q t3 

e i6/2 qu + e~ ie/2 qt4 



V 

0 

0 

e ^/2q- 3 + e - i S/2q t3 

e ie/2 ql 4 + e~ i9/2 q% l 

) 


(5.69) 


Because of Eq. (5.68), [Q(e 1 , £ 2 , v, 0 )] 2 is referred to as the amplification matrix of the a-e-4 scheme per 
marching cycle (i.e., per At) while Q(ei,e2,v,9) is referred to as the amplification factor per marching step 
(i.e., per At/2). Also by using Eq. (5.68), one has 

q*{n + m,9) = [Q(e 1 ,e 2 ,v,9)] 2m q*(n,0) m= 1,2,3...; n = 0, ±1/2, ±1, ± = 3/2, . . . (5.70) 

There is an important relation between the matrix Q(e 1 , 62 , v, 9) and the amplification matrix per marching 
step of the a-e scheme, i.e., 

M(e, v, 9) d =e ie/2 M_(e, v) + e~ i9/2 M + (e , v) = 

/cos(0/2) — ivsh\(9/2) — i(l — v 2 ) sin(0/2) \ — 00 < e < +00 (5.71) 

\ f(l — e) sin( 6 l / 2 ) (2e — 1) cos(0/2) — ivsin{9/2) J 
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In fact, it follow from Eqs. (5.63), (5.64), (5.69), and (5.71) that 


' e*/ 2 gn + e~ ie / 2 q+ e i0 / 2 qi 2 + e~ i0 / 2 q+ ' 
+ e- i0 / 2 q+ e i0 / 2 q 22 + e~ i0 / 2 q+ , 


= M(e2,v,6) 


(5.72) 


and 


e ^/2 q - 3 + e -w/2 g + e iei 2 g - + e -w/2 g + 
+ e- ie/2 9 4 + 3 e* /2 <& + e“ ie/ V 4 


j = M(c\,v, 9) 


(5.73) 


i.e., the 2x2 submatrices in the top-left and bottom-right corners of the matrix Q(ei, € 2 , is, 9) are M(e 2, 9) 

and M(ei,u,6), respectively. As such, the amplification matrix of the a-e scheme is identical to each of the 
above 2x2 submatrices of Q(e\, £ 2 , is, 6 ) except that the parameter e in the former matrix is replaced by £2 
or £1 in the latter two submatrices. By using this property, it will be shown that, for each set of £ 1 , £ 2 , is, 
and 9, the eigenvalues of Q(e\,e2,is,9) can be expressed as X±(ei,is,9) and A±(e 2 ,v,9) where the function 
A± ( e , v, 9) was defined in Eq. (2.21 ). 

To prove the above proposition, note that, for each set of £i, £ 2 , is, and 9, the eigenvalues A of 
Q(ei, £ 2 , is, 9) are the four roots of the characteristic equation 


det [<2 (£i, e 2 , is, 9) — XI4) = 0 (5-74) 

Hereafter /„ denote the n x n identity matrix. Let O 2 denote the 2x2 zero matrix. Then according to 
Eqs. (5.69), (5.72), and (5.73), the 2x2 submatrices in the top-left, bottom-left, and bottom-right corners 
of the 4x4 matrix [Q(ei, e 2 , is, 9) — XI4] are [M(e 2 , is, 9) — A/ 2 ], O 2 , and [M(ei,is, 9) — A/ 2 ], respectively. As 
such, by using elementary theory of determinant, we have 

det [<3(ei, e 2 , is, 9) - A/ 4 ] = det \M(e\,v,9) - X / 2 ] • det [M(e 2 ,i/,0) - A I 2 \ (5.75) 

It now follows from Eqs. (5.74) and (5.75) that, for each set of ei, £ 2 , and is, the eigenvalues of Q(ei,e 2, is, 9) 
are the four roots of the characteristic equations 

det [Mia, is, 9) - XI 2 ] = 0 (5.76) 

and 

det [Af (e 2) is, 9) - A/ 2 ] = 0 ((5.77) 

Because the roots of Eq. (5.76) and ((5.77) are A±(ei ,is,6) and A±(e 2 ,^, 6), respectively, the proof is com- 
pleted. QED. 

Because of the relation Eq. (5.70), for a given set of ei, e 2 , and is, the a-e- 4 scheme is said to be stable 
relative to this set if and only if, for all 9, —00 < 9 < + 00 , all the elements of the matrix [Q(ei,e2,is,9)] m 
remain bounded as m — > + 00 . By using the Jordan canonical form theorem [ 66 ] and the proposition we have 
just proved, it can be shown that, for each given set of ei, e 2 , is, and 9, all the elements of [Q(e\,e2,is,9)] m 
remain bounded as m — > +00 if and only if 


max {|A + (ei, is, 0 )| , |A_(ei,i',0)| , ]A+(£ 2 , is, 9)\ , |A_(e 2 ,i',0)|} <1 if Q(ti, e 2 , is, 9) is nondefective 


(5.78) 


and 


max {|A + (ei, v, 9) \ , \X-(ei,is,9 )\ , \X + (e 2 ,is,9)\ , |A_ (e 2 , is, 0)]} <1 if Q(ei, e 2 , is, 9) is defective (5.79) 

Note that an n x n matrix is said to be defective if and only if the number of its linearly independent 
eigenvectors is less than n [ 66 ] . 
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From the above discussions, it is seen that, for a given set of ei, £ 2 , and v, a necessary condition for the 
stability of a-e-4 scheme is 


max { | A_|_ (ei , v, 6)\ , |A_(ei,v,0)| , |A+(e 2 , v, 9)\ , |A_(e 2 ,i/,0)|} < 1 (5.80) 

— oo<c/<+oo 

By using Eq. (2.21) and a line of arguments which was used in establishing the stability condition Eq. (3.14) 
of [5] (Note: In [5], the a-e scheme is not defined for the case \u\ = 1), one can show that Eq. (5.80) 

0 < ei < 1, 0 < e 2 < 1 and \v\ < 1 (5.81) 

Thus Eq. (5.81) represents a set of necessary stability conditions for the a-e-4 scheme. 

Because Q(e i,e 2 ,y, 0) does become defective for some special values of ei e 2 , v and 9, according to 
Eq. (5.79), Eq. (5.80) and therefore Eq. (5.81) generally is not sufficient to insure the stability of the a-e-4 
scheme. However, it is expected that the necessary and sufficient stability conditions (not yet established) 
for the a-e-4 scheme would be formed by Eq. (5.81) and some exclusion conditions involving some boundaries 
of the stability domain defined in Eq. (5.81). 

From the above discussions, one concludes that the 2nd-order and 3rd-order terms which appear on the 
right sides ofEqs. (5.34) and (5.52) (which are associated with the 2x2 submatrices in the top-right corners 
of the matrices Q+ and Q-) play only a marginal role in the stability of the a-e-4 scheme even though they 
play a decisive role in the 4th-order accuracy of the scheme 

Next we will present a set of non-reflecting boundary conditions with 4th-order accuracy. 

5.5. Non-reflecting boundary conditions 

Let (j, n ) be a mesh point on the right spatial boundary. It has been shown numerically the current 
extension of Eq. (2.44) i.e., 

u] = u^Zx/li («*)" = [uxTjZyl, (uxx)j = (uxx)jl$, and (u SS x)] = [uxxxTjZyl (5-82) 

and a similar extension of Eq. (2.45) also form sets of non-reflecting boundary conditions. However, these 
simple extensions are not 4th order in accuracy. In the following, we present non-reflecting boundary 
conditions with 4th-order accuracy. 

Let (j, n) be a mesh point on the right spatial boundary. Again we assume 

(u*S*)? = (Uxxx^Zl/2 (5.83) 

However, («$$)" will be determined in terms of the mesh variables at (j — 1/2, n — 1/2) by imposing the 
conservation condition Eq. (5.13) which links the mesh variables at the two diagonally opposite mesh points 
(j, n) and (j — 1/2, n— 1/2). Because Eq. (5.13) <t=> 

(1 + v) [Uxx - (1 - v)u SS x] n j =(\ + v) [u m + (1 - ^UxxxZjZyl (5.84) 


we assume 


[u X x (1 ^ / )'U'xxx]j — [ 'U'xx T (1 v)u X xx\j_ij 2 (5.85) 

(Eq. (5.85) is a condition stronger than Eq. (5.84)). By combining Eqs. (5.83) and (5.85), one has 

( / Uxx')j ~ [Uxx “b 2(1 v)Uxxx\j_i/2 
Next («i)“ will be evaluated using the Taylor’s expansion 

-i n- 1/2 


(5.86) 


(«*)? = 


AX At 1 ( AX N 2 AX At 

-YUxx + Y u xt + 2 (-y) Uxxx + 


1 (At 


.. ^ xxt i i j u xtt 


(5.87) 


4 3- 1/2 
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With the aid of Eqs. (2.9), (5.3), and (5.17), Eq. (5.87) can be rewritten as 


(us) 


n 

j 


a x T 2(1 ^)wxx T 2(1 


^) a X xx\ 


n- 1/2 
3- 1/2 


(5.88) 


The mesh variable w" will be evaluated in terms of the mesh variables at (j — 1/2, n — 1/2) by imposing 
the conservation condition Eq. (2.8). Because Eq. (2.8) <t=> Eq. (5.30), we will assume Eq. (5.30a) (a stronger 
condition than Eq. (5.30)). As such one has 


n N 2(1 — v + v 2 ) 

(1 - v)Ux Uxx 


(l-z/)(l + zA) 


J 3 


+ {s-)7-vt 


(5.89) 


where is defined in Eq. (5.28). With the aid of Eqs. (5.83), (5.86), and (5.88), u" can now be 

evaluated in terms of the mesh variables at ( j — 1/2, n — 1/2). 

A similar treatment can also be applied at the left spatial boundary. As will be shown, these new 
non-reflecting boundary conditions are 4th-order in accuracy. 

6. The c-t-4 and c-r*-4 schemes 

The c-t-4 and c-t*-4 scheme will be built from the c-t scheme in a procedure similar to that was used 
to built the a-e-4 scheme from a-e scheme. Thus the basic equations including Eqs. (2.4), (2.9), (5.1)-(5.7), 
(5.11), and (5.17) will be assumed here. In addition, both the conservation conditions Eq. (2.13) and (5.18) 
will also be assumed here. 

6.1. The c-t-4 scheme 

Because Eq. (5.18) is equivalent to Eq. (5.14) assuming Eqs. (5.7), (5.11), and (5.17), in the current 
development, the mesh variable (u$s)" will be evaluated using Eq. (5.14). Moreover, because (i) Eq. (5.31) 
is equivalent to Eq. (2.13) assuming Eqs. (2.4), (2.9), (5.4), and (5.17); and (ii) Eq. (5.34) can be obtained 
from Eqs. (5.31) and (5.14), the mesh variable u " will be evaluated using Eq. (5.34). As such, the current 
schemes share with the a-e-4 scheme two component equations, i.e. , Eqs. (5.14) and (5.34). 

In the a-e-4 scheme, (u$x$)j is evaluated using Eq. (5.26) which is the substitution form of Eq. (2.19a) 
(by which (u 5 )" is evaluated in the a-e scheme) except that the symbol e in Eq. (2.19a) is replaced by ei 
in Eq. (5.26). Similarly, in the c-t-4 scheme, (u$s®)" is evaluated using the substitution form of Eq. (3.11) 
(by which (u®)" is evaluated in the c-t scheme) with the understanding that the symbol t in Eq. (3.11) is 
replaced by ti in the substitution form. Thus we have 


(u X xx 


1 

2(1 + Tl) 



(1 + ‘Iv 


7~l)U'XXX 


,n-l/2 

b+1/2 


\a X x T (1 2i/ 



(j,n) £ fi; ti ^ - 1 


( 6 . 1 ) 


In the c-t-4 scheme, (u$)" will be evaluated using an equation which is the result of the following 
assumptions: (i) (u$)" = [(4 (t 2)]" (see Eq. (5.48)); ( ii) the mesh variable ( u^xx )” which appears on the 
right side of Eq. (5.48) is evaluated using Eq. (6.1); and (iii) the parameter t which appears on the right 
sides of Eqs. (5.44) and (5.45) is replaced by T 2 . Then, by using Eqs. (5.44), (5.46), (5.48), and (6.1), one 
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has 


(Ux) 


n 

3 


1 

2 


u 

1 + t 2 


1 - t 2 + 2i/ 1 r(l — t 2 + 2i/) 2 

Up- ~f“ — 

1 + r 2 x 2 L l + r 2 


(l + r 2 ) 2 l 

3(1 + Tl)- 


'Uxx 


1 r (l + t 2 ) 2 (1 - ti + 2i/) 

6 1- 1 + Tl 


(1 — r 2 + 2t/) 3 ' 
1 + t 2 


Uxxx 


■x n— 1/2 

J j+1/2 


1 

2 


u 

1 + r 2 


1 — r 2 — 2i/ 1 r (1 — r 2 — 2v) 2 

1 + r 2 * 2 L 1 + r 2 


(l+r 2 ) 2 ] 

3(1 + Tl) . 


^ XX 


( 6 . 2 ) 


1 r(l + t 2 ) 2 (1 - ti - 2v) 


6 L 


1 + Tl 


(1 - t 2 - 2vf - 
1 + t 2 


Uxxx 


■x n- 1/2 
J J — 1/2 


(j,n) € f2; Tl ^ -1; t 2 ^ -1 

By definition, the c-t-4 scheme is formed by Eqs. (5.14), (5.34), (6.1), and (6.2). It enforces the local 
conservation relations Eqs. (2.13) and (5.18). According to an argument given in Comment (d) presented at 
the end of Sec. 2.1, it will also enforce a global conservation relation for each of the fluxes of h* and h* x over 
any space-time region that is the union of any combination of CE(j, n), (j, n) £ Cl. Moreover, the scheme is 
completely explicit, i.e., the four independent mesh variables at each ( j , n ) can be evaluated directly in terms 
of those at the (n — l/2)th time level without resorting to solving a system of linear equations involving the 
independent mesh variables at ( j,n ) such as Eqs. (2.7a) and (2.8a). 

The c-t and a schemes become the same scheme when v = t = 0 (see Eq. (3.12)). It happens that 
the c-t-4 and a-4 schemes also become the same scheme when v = ti = t 2 = 0. In fact, for the case 
v = Ti = t 2 = 0, Eqs. (6.1) and (6.2) reduce to 


(n X xx 


)" 

h 



XX ^ 


XXX 


xn- 1/2 
>3 + 1/2 


[n X x 


+ u 


XXX 


xn-1/2 
> 3 - 1/2 


(6.3) 


and 


1 

v 1 7 11-1/2 

( 1 N 

. 77—1/2" 

o 

1 u-u x + -u xx I 

[u + u x + -u xx 


Z 

V 6 J j+ 1/2 

\ 3 y 

3~ 1/2. 


(6.4) 


respectively. It can be shown that Eqs. (5.26) and (5.52) (the counterparts of Eqs. (6.1) and (6.2) in the 
a-e-4 scheme) also reduce to Eqs. (6.3) and (6.4), respectively, when u = ei = e 2 = 0. Because (i) the a-4 
scheme is a special case of a-e-4 scheme with ei = e 2 = 0 and (ii) Eqs. (5.14) and (5.34) are shared by the 
c-t-4 and a-e-4 schemes, the proof is completed. QED. 

Because the conservation relations Eqs. (2.7), (2.8), (5.12), and (5.13) are satisfied by the a-4 scheme for 
all v , it now follows that these conservation relations are satisfied by the c-t-4 scheme when v = ti = t 2 =0. 

The c-t and c schemes (the special case of the a-e scheme with e = 1/2) become the same scheme when 
t = 1 (see Eq. 3.12). Let the c-4 scheme be the special case of the a-e-4 scheme with ei = e 2 = 1/2. Then 
one can also show that the c-t-4 and c-4 schemes become the same scheme when Ti = t 2 = 1. In fact, for 
the case n = t 2 = 1, Eqs. (6.1) and (6.2) reduce to 


( v ---) n — I 17,, 2nu \ n ~ 1 / 2 _ ( v 

\ u, xxxjj — ^ y^xx Uj xxx)j_ |-i/2 yu'xx ^xxx J j—i/2 


(6.5) 


and 


(«i)" = Y2 { [ 3u “ 6l/u x + 2(3^ 2 - 1 )u xx + 4 i/(l - v 2 )u xxx ]™ + l j* 
- [3u -6 vu x + 2(3 v 2 - 1 )u xx + 4 i/(l - v 2 )u xxx \ 


( 6 . 6 ) 
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respectively. It can be shown that Eqs. (5.26) and (5.52) also reduce to Eqs. (6.5) and (6.6), respectively, 
when ei = £2 = 1/2 (i.e., when the a-e-4 scheme becomes the c-4 scheme). QED. 

6.2. Matrix form of the c-r-4 scheme 

The component equations of the c-r-4 scheme, i.e., Eqs. (5.14), (5.34), (6.1) and (6.2), can be cast into 
the matrix equation 

q{j,n) = Q-(ti,T 2 , v) q(j + 1/2, n- 1/2) + Q + (n , r 2 , v) q{ j - 1/2, n- 1/2) (, j,n ) G Q (6.7) 

where q(j,n) is defined in Eq. (5.54), and 


with 


and 
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def 
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-oi 
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V 0 
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qf A / 
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(l + r 2 ) 2 

7T' 

923 

== dz 

12(1 + n) 

(1 

- r 2 

+ 2 i /) 3 



12(1 + 75) 




6 


-F 


(1 - t 2 + 2+) 2 
4(1 + r 2 ) 


934 = ± 


") 943 — 


def + 1 


2(1 + Ti) 


i 944 — 


def (1 T\ + 2v) 


2(1 + Ti) 


(6.8) 


(6.9) 


(6.10) 


( 6 . 11 ) 


To further explore the relations between the c-r-4 and c-r schemes, note that the component equations 
of the c-r scheme, i.e., Eqs. (2.10) and (3.11), can be cast into the matrix form 

g(j, n) = M_(r, v) q{j + 1/2, n - 1/2) + M+(r, i/) g(j - 1/2, n - 1/2) (j, n) G ft (6.12) 

where q(j,n) is defined in Eq. (5.61), and 


M±(t,v) d = i 


( l±v ±{l-v 2 ) 
+1 1 - r 3p2v 

V 1 + r 1 + r 


With the aid of Eqs. (6.9)-(6.11) and (6.13), it can be shown that 

= M ± {t 2 ,v) 
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9n 

9l2 

Z~|- 


921 

922 



933 

934 


-oi 

943 

944 


= M±(n,u) 


— 00 < T < +00 


(6.13) 


(6.14) 


(6.15) 
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By comparing Eqs. (6.14) and (6.15) with Eq. (6.8), one now arrives at the conclusion that the 2x2 
submatrices in the top-left and bottom right corners of the 4x4 matrix Q±(ti 5 T 2 , is) are M±(t 2 , is) and 
respectively. 

6.3. Stability conditions of the c-r-4 scheme 

By using similar arguments presented in Sec. 5.4, one concludes that the amplification matrices per 
marching step of the c-r-4 and c-r schemes are 


Qin,T2,is,9) = e l9/2 Q-iT 1 ,T 2 ,is) + e l9/2 Q+iTi 

,T2,1S) 




/ e i9 / 2 q 41 + e~ w / 2 q+ e i9 / 2 q 42 + e~ t9 / 2 q+ 2 

e ie/2^ + e -ie/2^ 

e i9/2 q 14 + e i9/2 qt 4 

\ 


e ie ' 2 q 21 + e~ w / 2 q+ e i9 ' 2 q 22 + e~ i9 / 2 q+ 2 

e iei2^ + e -ioi2^+ 

e i9 ' 2 q 24 + e~ i9 ^ 2 q 24 



0 0 

e i9 / 2 q 33 + e~ i9 / 2 q+ 

e i9l2 q 34 + e~ i9 / 2 qt 4 



\0 0 

e i9 ' 2 q 43 + e-'°/ 2 qt 3 

e i9 / 2 q 44 + e- iS / 2 q+ 4 

/ 


and 

M(r, is, 9) d =e ie / 2 M_(T, is) + e~ l9 ^ 2 M + ij, is) = 


(6.16) 


/ cos(0/2) — ivsm{6/2) —i( 1 — zA)sin(#/2) 

tsin(0/2) 

V 1 + r 

respectively. Moreover, by using Eqs. (6.14)-(6.17), one concludes that 

j = M{t 2 ,is,6) 


e ie /% t + e~ i9 / 2 q+ e ie / 2 q 42 + e~ i6 / 2 q+ 


, e"/ 2 ^ + e-«/*q+ e i9 / 2 q 22 + e~ i9 / 2 q+ 


(1 — r) cos(0/2) + 2ii/sm(0/2) 
1 + r 


(6.17) 


(6.18) 


and 


e i9 ' 2 q 33 + e~ i9 / 2 q+ e i9 /% 4 + e"*/ 2 g+ 

e i9/2 q 43 + e~^/ 2 q+ e iS / 2 q 44 + e~ i9 / 2 qt 4 


j =M(t u v,9) 


(6.19) 


i.e., the 2x2 submatrices in the top-left and bottom-right corners of the matrix Q(t\ , T2, is, 9) are M(t2, is, 9) 
and M(ti,u,9), respectively. As such, the amplification matrix of the c-t scheme is identical to each of the 
above 2x2 submatrices of Q(t 4 , T 2 , is, 9) except that the parameter t in the former matrix is replaced by T 2 
and Ti in the latter two submatrices. By using this property and using a line of arguments similar to that 
presented in Sec. 5.4, one can show that, for each set of n, r 2 , is, and 9, the eigenvalues of Q(n, r 2 , is, 9) are 
the union of those of M(t\, is, 9) and M(j 2 , is, 9). As a result, by using a line of arguments presented in the 
last part of Sec. 5.4, one can infer from the stability conditions Eqs. (3.13) of the c-r scheme that 


is 2 <1, n > Toils 2 ), r 2 > T 0 {is 2 ), 


i^ 2 , Ti) yf (1, 1), and iis 2 , r 2 ) ^ (1, 1) 


( 6 . 20 ) 


form a set of necessary stability conditions for the c-r-4 scheme. 

6.4. The c-r*-4 schemes 

Let hiis) and h 2 is) be any two functions that meet the same requirements imposed on the function ft.(s) 
introduced in Sec. 3.2. Then, because of the close relations between the c-r and c-r-4 schemes established 
above, here a c-r*-4 scheme is defined to be a special c-r-4 scheme with 

Ti = h\iis 2 ) and r 2 = h 2 iis 2 ) i\is\ < 1) (6-21) 
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According to Eq. (3.22), one can assume h\(s) = fi 2 (s) = \fs- As such a special c-r*-4 scheme can be defined 
by the relation 

ti=t 2 = H (M <1) (6.22) 

For robust real-world applications, we may also introduce a special c-r-4 scheme which is defined by an 
extension of Eq. (3.24), i.e., 

n = Pi \v\ and r 2 = /3 2 M (Pi > 1; Pi > 1; M < 1) (6.23) 

where /3i > 1 and p 2 > 1 are two adjustable parameters. 

7. Fourth-order solvers for the inviscid Burgers equation 

Consider the inviscid Burgers equation 


du df(u) 
dt dx 

(7.1) 

where 


\ def u 2 

f(u) = y 

(7.2) 

Because, in E 2 , Eq. (7.1) can be expressed as V • h = 0 with 


h = (f,u) 

(7.3) 


Gauss’ divergence theorem in the space-time E 2 implies that Eq. (7.1) is the differential form of the integral 
conservation law Eq. (2.2) if h is defined by Eq. (7.3). Moreover, for a smooth function u(x,t), Eq. (2.2) 
implies the conservation condition Eq. (5.9) if 


r def = ( cPf_ d 2 u \ 

xx ~ dx 2 ~ \dx 2 ’ dx 2 J 


(7.4) 


7.1. The inviscid Burgers c-r-4 scheme 

In this section, the c-r-4 scheme will be extended to become a solver for the inviscid Burgers equa- 
tion. Here again Eq. (5.1) is assumed. However, for any (x,t) £ SE (j,n), the flux function f(u) will be 
approximated by 

f*(x,t,j,n) d = /j 1 + (f x )j(x - Xj) + (ft)j(t - t n ) + ^(f xx )](x - Xj ) 2 

+ (f xt )](x - Xj)(t - n + - n 2 + ^ - Xj ) 3 (x,t) £ SE (j,n) (7.5) 

+ \(u x t)](x Xjf(t n + ±(f xtt )](x - Xj)(t - n 2 + l(fut)](t - n 3 


Here /?, (f x )J, (ft)], (fxx)], ( fxt )", ( ftt )" (fxxx)], (fxxt)], ( fxtt )], and ( fttt )" are constants in SE(j,n). 
They can be considered as the numerical analogues of the values of /, df/dx, df /dt, d 2 f/dx 2 , d 2 f/dxdt, 
d 2 f/dt 2 , d 3 f/d x 3 , d 3 f/dx 2 dt, d 3 f/dxdt 2 , and d 3 f/dt 3 at the mesh point (j, n) , respectively. Thus 
f*(x,t,j,n) represents a 3rd-order Taylor’s approximation of /. Moreover, as will be shown, by impos- 
ing proper assumptions, the coefficients on the right sides of Eqs. (5.1) and (7.5) can be explicity expressed 
as the functions of the independent unknowns u ”, (u x )], (u xx )] and (u xxx )]. 

For any (j,n) £ O, let u = u*(x,t; j,n) and f(u) = f*(x,t,j,n) satisfy Eq. (7.1) uniformly in SE (j,n), 


i.e., 


du*(x, t ; j, n) df*(x,t,j,n) = 
dt dx 


(j, n) £ n 


(7.6) 
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It can be shown that Eq. (7.6) <t=>, for all (j,ri) £ f 2, 

(«*)” = 

( )] = ~(/xx)? 

(«**)” = -(/**)" 

(Uxxt)" = 

0 «xtt )] = -(/xxt)" 

and 

(wtit)" = — (/lit)™ 

Moreover, by taking the numerical analogues of all the differential forms 
one has, for all (j, n) £ fi, 

y? = |(« 2 )? 

(/*)? = (««*)? 

(ft)] = (uu t )] 

(fxx)j = [(uu xx + (Mo;) 2 ] " 

(fxt)j = ('U'U'xt “t” UxUt) j 

( fu )]=[(uu tt + (« t ) 2 ]; 

(fxxx)j = (’U'Uxxx + 3 U X U XX ) j 
( fxxt)j = (Wxxt “t” 'U’t’Uxx “t” 2u x U x t.) j 
(fxtt)j = ('UUxtt “t” ^x^tt “t” 2 UtUxt) j 

( fttt. )] = ( uuttt + 3 ut.utt)] 

By substituting Eqs. (7.14), (7.16), and (7.19) into Eqs. (7.7), (7.8), 

(«*)? = -K)" 


(7.12) 

of Eq. (7.2) up to the third order, 


(7.19) 

(7.20) 

(7.21) 

(7.22) 

and (7.10), respectively, one has 


('U'xt) j — 'U’U X x 4“ (^-x) 


(^ xxt) j — 1 3'Ux'U'xx) j 

respectively. In turn, by substituting Eqs. (7.23)-(7.25) into Eqs. (7.15), (7.17), and (7.20), one has 

(ft)] = -( u 2 u x )] 

(fxt)j = - [2u(u x ) 2 + U 2 ^] " 

and 

(fxxt)j = [2(ria,) + Quu x u xx u u xxx ^ j 

Moreover, by substituting Eqs. (7.27) and (7.28) into Eqs. (7.9) and (7.11), respectively, one has 

= [2m(m£c) 2 + w 2 Mxx]” 
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and 


(7.30) 


{^xtt ) j — [2(w^) Qu'Ux'U'xx "1” U U X xx\ j 

respectively. 

Next, by substituting Eqs. (7.23), (7.24), (7.29), and (7.30) into Eqs. (7.18) and (7.21), one obtains 

( fa )j = [3 u 2 {u x ) 2 +v?u xx ] n . (7.31) 

and 

(fxtt)j — [0'u(u x ) T- 9 u u x iixx 4- u u xxx (7.32) 
In turn, by substituting Eq. (7.32) into Eq. (7.12), one has 

(^ ttt)j = [0'u('u a; ) t 9 u u x u X x 4“ u u X xx\j (7.33) 

Finally, by substituting Eqs. (7.23), (7.29), and (7.33) into Eq. (7.22), one has 

{fttt )j = - [12 u 2 (u x ) 3 + 12 u 3 u x u xx + u 4 u X xx ]” (7.34) 

According to the set of sixteen equations formed by Eqs. (7.13), (7.14), (7.16), (7.19), and (7.23)- 
(7.34), all the twenty coefficients on the right sides of Eqs. (5.1) and (7.5) other than u”, ( u x )", ( u xx )", and 
( UxxxYj are explicit functions of the latter coefficients. Thus there are only four independent mesh variables 
associated with each (j, n) € O. In fact, by substituting these sixteen equations into Eqs. (5.1) and (7.5), 
one has 


“(x,t;j,n ) =Uj + (u x )]{l - - t n ) + [(u x )™] 2 (t - t n ) 2 j [a; - xj - u™(t - t n ) 


( a xx) 


(^XXx) i 


(7.35) 


+ [i - 3 (« x )"(t - n] [x - xj - u](t - n] “ + v : [® - Xj - u]{t - 1 


and 


n\2 


f*(x, t,j, n) =^~y~ + ("x )'/{"'/ + [a: - xj - 3 u](t - t n )\ 

- [{u x )j] 2 (t - t n )[x - Xj - 2 Uj(t - t n )] | [x - Xj - Uj(t - t n )] 
+ ( "i[ )f + («x)i [x - Xj - 4u?(t - t n )\ } [a; - xj - u](t - t n )] 2 


(u XX x) j 

— 


X — Xj — u ” ( t — t r ‘ 


In turn, by substituting Eqs. (7.35) and (7.36) into Eq. (5.5) and 

,* +. ■ \ d _£ f d2 f*(x,t,j,n) 

f xx \Xj t, j, n) fix' 2 


one has 


c(a ; ?t;j,u) — {Uxx) j T i^'xxx) j (x x j ) ('uu xxx 4“ 3 u x u xx )j{t t ) 


(7.36) 


(7.37) 

(7.38) 


and 

fx X (.Xit\j-)ri) = \uU X x 4” (^x) }j 4” (uU xxx 4~ 3u x U xx ) j (x Xj) [2(Ua,) 4“ ()UU x U xx 4” a (t t ) (7.39) 
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Note that, for the same reason presented immediately following Eq. (2.6), except for Bux evaluation, 
Eqs. (7.35) and (7.38) are applicable even for (x,t) £ SE(j,n). 

Next let the numerical analogues of Eqs. (7.3) and (7.4) be 

— * 

h*(x,t;j,n) = (f*(x,t,j,n),u*(x,t;j,n)) (x,t) G SE (j,n) (7.40) 

and 

Kx( x ’ t 'J’ n ) = (fL( x ’ t 'J’ n )> u *xx( x i t 'J, n )) (x,t) € SE (j,n) (7.41) 

respectively. Also, let 

(u*)? d = ^(O", (|) 2 (%x)?, and (u^)] = (u xxx )] (i,n)e Q (7.42) 


and 


def O At 


AX 


(«/*)? = 


def ( u x ) ; At 


AX 


(***)? = 


def ( u xx) 7 At 


AX 


and (vxxx)l = 


def ( u m)j At 


AX 


(j) n) & ft (7.43) 


Obviously, («$)", («sx)", and (u S xx)(j defined here are identical to those defined in Eqs. (2.9) and (5.17). 
However, the “Courant number” v™ and its “derivatives” (z's)™, (%)", and (v-cxx)™ are mesh-point depen- 
dent in the current development. 

Like the c-r-4 scheme, the current extension will be construct to enforce the least stringent conservation 
conditions Eqs. (2.13) and (5.18). In fact, by using Eqs. (7.38), (7.39), (7.41), (7.42), and (7.43), it can be 
shown that Eq. (5.18) <+> 


i- 1/2 
0 + 1/2 


[n X x) j — 2 { \ V xQ V x l)llx "t“ [1 ^(1 6v x )\Uxx (1 V )lla:xx] , 

A [+x(l 2 v^)Ux 4- [1 4~ ^(1 6r , x)]ifxx 4” (1 r' )wxxx] j_\/2 J* 

Moreover, by using Eqs. (7.35), (7.36), (7.40), (7.42), and (7.43), it can be shown that Eq. (2.13) <t=> 


(7.44) 


= {(l - 7j)u - [l - u 2 + 2vxV 2 - 4 (vx) 2 v 2 ] u x + [ 
+ (l + 2') u 4” [l — + 2 + 2 I'xi' 2 — 4 (r , $) 2 ^ _ ] + [ 


2(1 — v 3 ) 31 

Q H - 4z/ x ^ 1 

O LL XXX 

, n- 1/2 

3 J 

3 

J j + 1/2 

2(1 + ,3) 

4^1/ 1 U X x H - 

o 11 XXX 1 

n-1/2 

3 J 

3 J 

J — 1/2 


(7.45) 


Expressing the mesh variable (u^g)” on the left side of Eq. (7.45) in terms of the mesh variables at the 
(n — l/2)th time level by using Eq. (7.44), one obtains 


(2-1/) [3(1 — v 2 ) + 2(3jA — 1)(1 — 2vx)vx\ v(l - v 2 )(l - bux) (1 - v 2 ) 2 ^-1/2 

, U r . Mx Uxx H - Uxxx f 

4 6 3 6 ) j+ 1/2 


( 2 + is ) [3(1 - v 2 ) + 2(3iA - 1)(1 - 2ux)Vx\ v(l - i/ 2 )(l - bvx) (1 - v 2 ) 


~ 'U'xxx f 

6 Jj - 1/2 


-1/2 

-1/2 

(7.46) 
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It has been shown that, by imposing the conservation conditions Eqs. (2.13) and (5.18), one can obtain 
Eqs. (7.44) and (7.46) by which both u" and (uxx)™ can be explicitly evaluated in terms of the mesh 
variables at the (n — l/2)th time level. In the following, by using techniques developed earlier, similar 
explicit expressions will be developed for (u $ )” and (u^xx)"- 

To proceed, let P + and P~ be the points depicted in Fig. 3, except that the parameter r shown in the 
figure be replaced by the mesh-point dependent parameter (n)”. Then the x-coordinates of P + and P~ are 


x(P+) = Xj + [1 + (ri)”]— = x j+1/2 - [1 - (ri)”]— , -oo < (ri)^ < +oo (7.47) 

and 

At At 

x(P-) = Xj - [1 + (n)"] — = Xj _ 1/2 + [1 - (ti)"] — , -00 < in)] < +oo (7.48) 

respectively. In the current development, we assume 



! Uxx(P + ) Uxx(P ) 

[1 + (n)"] ax/2 

(7.49) 

where 

u xx (P+) ^ u* xx (x(P + ),t n ;j + 1/2, n- 1/2) 

(7.50) 

and 

Uxx(P-) = f u* xx (x(P~), t n ; j - 1/2, n - 1/2) 

(7.51) 


Because |P _ P+| = [l + (ti)™] az/2, the expression on the right side of Eq. (7.49) represents a central- 
difference approximation of (u xxx )'j. With the aid of Eqs. (7.38), (7.42), (7.43), (7.47), (7.48), (7.50), and 
(7.51), it can be shown that Eq. (7.49) 


id^xxx) i — 


.— . 1 / 2 " “ 71 

[(1 bv x )u x x (1 Ti + 2z ^U xxx (1 bh , x )u x x + (1 TL 2r')tixxx j 


- 1/2 
i- 1/2 


2[1 + (ti)3 


(7.52) 


To simplify notation, in Eq. (7.52) and hereafter in this section we adopt a convention that can be explained 
using two expression on the right side of Eq. (7.52) as examples, i.e. , 


[(1 bVx^Xlxx 


(1 T\ + ‘2i / ')u xxx 


n-l/2 

3+1/2 


[(1 - 6Vx)Uxz]"d/2 


[l-(r I )," + 2^-'/ 2 2 ](« 


\n-l/2 
>3 + 1/2 


and 


[(1 bVx^Uxx + (1 


7~1 — 2 v)u 


XXX 


n-l/2 

3 - 1/2 


[(1 - bv x )u, 


n-l/2 

xx I / — 1/2 


+ [l - (ri)" - (Uxxx 


\n— 1/2 
> 3 - 1/2 


In other words, the parameter ti which appears on the right side of an equation such as Eq. (7.52) 
and the mesh variable which appears on the left side are associated with the same mesh point 

(/»• 

Next, let P + and P be the points depicted in Fig. 3, except that the parameter r shown in the figure 
be replaced by the mesh-point dependent parameter ( 72 )”. Then the x-coordinates of P + and P~ are 

x(P+) = Xj + [1 + (r 2 )"]^ = x j+1/ 2 — [1 - (t 2 )"]^, -00 < (r 2 )" < +00 (7.53) 

and 

At At 

x (P~) = x 3 ~[ 1 + ( T 2 )]]—= x j- 1/2 + [1 - (t 2 )”]— , -00 < (r 2 )" < +00 (7.54) 
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respectively. In the current development, (u x )” is approximated with the current version of the expression 
on the right side of Eq. (5.46), i.e., 


, « u d *(P+)-u 4 *(P~) 

[x)j ~ [l + (r 2 )”] ax/2 


{[1 + (r 2 )"]Aa:} 2 
96 


(rr X xx 


(7.55) 


where 

u 4 *(P±) d £ u * {x{P ± ),t n \j ± 1/2, n — 1/2) (7.56) 

According to Eq. (7.35), u 4 *(P + ) is the 3rcl-order Taylor’s approximation (i.e., the error is O [(ax) 4 ]) of u 
at P + evaluated using the mesh variables at (j + 1/2, n — 1/2) while u A *{P~) is that at P~ evaluated using 
the mesh variables at ( j — 1/2, n — 1/2). As such, Eq. (5.43) implies that Eq. (7.55) is an approximation of 
3rcl-order accuracy for the value of du/dx at the mesh point ( j,n ). 

Next, with the aid Eqs. (7.35), (7.42), (7.43), (7.53), and (7.54), it can be shown that 


4* 


(■ P ± ) ={« + (1 - t 2 ± 2v)[l - 2u s +4 (Px) 2 }ux + ^(1 - r 2 ± 2 j /) 2 (1 - 6 i/ s )u s 

1 'l r 

~F g(l ^2 i 2z/) UxxxJ 


in- 1/2 

*j± 1/2 


(7.57) 


In turn, by using Eqs. (7. 42), (7. 52) and (7.57), it can be shown that Eq. (7.55) <t=> 


(us) 


n 

j 


2(1 + r 2 ) 


(1 - t 2 + 2v) [l - 2vx + 4(j/j) 2 ] Ux 

2(1 + r 2 ) 


+ (1 - 6ux) 


■(1 — r 2 + 2v) 2 

- 4(1 + T 2 ) 


(1 + t 2 ) 2 1 

12(1 + n)J 


'U'xx 


~ (1 + ^~2 ) 2 ( 1 ~ Ti + 2v) 
12(1 + n) 


(1 — r 2 + 2i/) 3 i 

12(1 + r 2 ) . 


rt X xx 


s. n— 1/2 

J j+1/2 


2(1 + r 2 ) 


(1 - r 2 - 2i/) [l - 2i>x + 4 (^ x ) 2 ]u x 

2(1 + t 2 ) 


+ (1 — 6%.) 


~ (1 - t 2 - 2v) 2 
- 4(1 + r 2 ) 


(l + ^2) 2 

12(1 + n) 


'U'xx 


- (1 + x 2 ) 2 (1 - Ti - 2^) 

12(1 + n) 


(1 — r 2 — 2i/) 3 - 

12(1 + r 2 ) . 


Hxxx 


n— 1/2 

J j-1/2 


(7.58) 

Note that the simplifying convention established earlier for the parameter t \ in Eq. (7.52) 
applies here for both n and r 2 , i.e., everywhere on the right sides of Eqs. (7.57) and (7.58), 
the symbols n and r 2 are the shorthands for (n)" and (r 2 )”, respectively. 

By definition, the Burgers c-r-4 scheme is formed by Eqs. (7.44), (7.46), (7.52), and (7.58). It enforces 
the local conservation relations Eqs. (2.13) and (5.18). According to an argument given in Comment (d) 
presented at the end of Sec. 2.1, it will also enforce a global conservation relation for each of the fluxes 
of h* and h* x over any space-time region that is the union of any combination of CE (j,n), ( j,n ) € fl. 
Moreover, even though it is a solver for a nonlinear equation Eq. (7.1), by explicitly specifying (n)” and 
(r 2 )" as functions of the mesh variables at ( j — 1/2, n — 1/2) and ( j + 1/2, n — 1/2), the scheme will become 
completely explicit, i.e., the four independent mesh variables at each (j, n) can be evaluated in terms of 
those at the (n — l/2)th time level without resorting to solving a system of nonlinear equations involving 
the independent mesh variables at (j. n). Also note that, because v/ = (z/ s )” = 0 when At = 0, Eqs. (7.52) 
and (7.58) reduce to Eqs. (6.3) and (6.4) respectively when At = n = r 2 = 0 In fact, it can be shown that 
the Burgers c-r-4 scheme satisfies Eqs. (2.7), (2.8), (5.12), and (5.13) (i.e., the conservation relations over 
each CE±(j, n)) when At = n = r 2 = 0. 
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To turn the current scheme into an explicit scheme, as an example, we can assume that 

(n)" = Pi\v™\ and ( 75 )" = feWj I (/?i > 1; fe > 1) (7.59) 

where (3\ > 1 and /3 2 > 1 are two adjustable constant parameters. According to Eqs. (7.43) and (7.46), 
Eq. (7.59) effectively defines (n)” and (r 2 )" in terms of the mesh variables at (j + 1/2, n — 1/2) and 
(j — 1/2, n — 1/2). As another example, one can assume 

(n)" = /?i max { | I, ^"+i/ / 2 2 |} ancl ( T 2 )" = fc max { | i^Zl /2 I, l^+i/^l} (ft > !; ^2 > 1) (7.60) 

where /3i > 1 and /3 2 > 1 are again two adjustable constant parameters. Other possible choices of (ri)" and 
(r 2 )" are specified in Eqs. (7.108)— (7.110). 

In the following, it will be shown that the stability conditions of the Burgers c-r-4 scheme are essentially 
similar to those of the c-r-4 scheme. 

7.2. Stability conditions of the inviscid Burgers c-r-4 scheme 

As explained in Sec. 5.4, for linear solvers such as the a-e-4 or c-r-4 scheme, propagation of round- 
off errors is governed by the same linear system of equations governing the associated independent mesh 
variables. As will be shown, the above simple relation is not valid for a nonlinear solver such as the Burgers 
c-r-4 scheme. 

To proceed, note that Eqs. (7.1) and (7.2) imply that the dependent variable u has the dimension of 
At /ax. Thus i/", (i^x)", (%)", and {vxxx)/j defined in Eq. (7.43) are all dimensionless. As a result, by 
multiplying both sides of each of Eqs. (7.44), (7.46), (7.52), and (7.58) with the factor At/ Ax, one can turn 
the latter equations into their dimensionless versions. The dimensionless versions can be obtained from 
the original versions by replacing the symbols it, Ux& and Uxxx with v, Vxxi and v xxx-> respectively. 
Explicitly, the dimensionless versions of Eqs. (7.46), (7.58), (7.44), and (7.52) are 


„ f(2 -v) [3(1 — u 2 ) + 2(3zA — 1)(1 — 2ux) v x\ z/(l - n 2 )(l - 6i/ s ) 

^ l 4 V 6 ^ + 3 

f(2 + i/) [3(1 — v 2 ) + 2(3zA — 1)(1 — 2t/x)i>x] v{\ - v 2 ){l - 6u s ) 

+ l 4 6 3 


(l^ 2 ) 2 

6 

(1~^ 2 ) 2 

6 


^ n- 1/2 
> 3+1/2 
^ n- 1/2 

J J — 1/2 


(7.61) 


(**)? = 


2(1 + r 2 ) 


(1 - r 2 + 2 v) [l - 2vx + 4(ti $ ) 2 ] Vx 
2(1 + r 2 ) 


+ (1 - 6^x) 


(1-t 2 + 2v) 2 (1 + r 2 ) 2 


4(1 + 75) 12(1 + n). 


' (1 + r 2 ) 2 (l - ri + 2v) 
12(1 + n) 


(1 — r 2 + 2n) 3 i 
12(1 + 75 ) - 


Uxxx 


n- 1/2 

J j+1/2 


2(1 + r 2 ) 


(1 - r 2 


2v)[l - 2vx + ^{vx) 2 ]^x c r(l — r 2 — 2v) 2 

2(1 + r 2 ) + (1 - 6, '* ) 4(1 + t 2 ) 


(l + ^2) 2 
12(1 + n) 




~ (1 + t 2 ) 2 (1 - n - 2v) 
12(1 + n) 


(1 - r 2 - 2n) 3 ' 
12(1 + 75 ) - 


Uxxx 


A n-1/2 

J J-l/2 


(rf 



l)(i'$) 2 + [1 - K 1 - 6z/ s )]z/ ss 


(1 ) Uxxx 


n-1/2 
3+ 1/2 


+ [(1 — 2jl 5 ;.)(ti , $ ) 2 + [1 + ^(1 — 6zi'2)]j y ££ + (1 — U 2 )u S xx j_ 1/2 | 


(7.62) 

(7.63) 
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and 


(u—' l 1 ? — 
y^xxxjj — 


(1 - 6Vx)Vxx - (1 - n + 2v)v S xx_ 


1 ”+l/2 “ [ (1 - 


T (1 T\ ‘2h / ')V X xx 


n—1/2 
3~ 1/2 


2[l + (n)y] 


(7.64) 


respectively. 

In the current study, we assume that 


1 ( Uxxx )? | « | («**) ] 1 « 1 («*)? | « | U] | (j 

i, n ) G H 

(7.65) 

Wj\ < 1 (j, n) G 0 


(7.66) 

and 



l(n)"| = O(0”|) and |(r 2 )?| = 0(\v?\), 

0, n) G O 

(7.67) 

Note that, by using Eqs. (7.35) and (7.42), one has 



At II 

u*( Xj ± —,t n ;j,n) = u] ± (u s )] + ~(uxx)3 ± g(te)? 

0, n) G ft 

(7.68) 


Thus, Eq. (7.65) simply states that the magnitude of a term on the right side of Eq. (7.68) is much larger 
than that of the next higher-order term to its right. In view of Eqs. (7.42), obviously, this assumption can 
be made valid by choosing a small enough value of ax. Moreover, by using Eqs. (7.43), it can be shown that 
Eqs. (7.65) and (7.66) ++ 


|(i'***)"| <| (v*s)"| <| (vx)j\ <K"I < 1 G n (7.69) 

With the aid of Eqs. (7.42) and (7.43), it is seen that the assumption |(^$)"| <C| i/"| in Eq. (7.69) <t=> 

Ax(u x )j < Uj ( j,n ) G (7.70) 


Within a smooth solution region , (tt x )" is a good approximation of the value of du/dx at ( j,n ). Thus 
Eq. (7.70) implies that 

l u "+i /2 - Vi/2 2 1 “"± 1 / 2 1 0> n) e O (7.71) 

By using Eq. (7.43), one can cast Eq. (7.71) into its dimensionless form, i.e., 


I n—1/2 n— 1/2 1 i n— 1/2 1 

g j+1/2 V j-1/2 I ^ I U j±l/2 I 


O', ™) g n 


In turn, Eq. (7.72) implies that 


lO? - *"i « W 


0, n) G O 


with 

Eq. (7.73) implies that 


def 


1 ( n-1/2 , n-l/2\ 

2 0+1/2 + *0-1/2 J 


0, n) G n 


n-l/2 

*0+1/2 *0 


0, n) G O 


(7.72) 


(7.73) 

(7.74) 

(7.75) 


Hereafter the symbol “=” denote the fact that the magnitude of the difference of the quantities on the both 
sides of the symbol is much less than the magnitude of either quantity itself. In other words, either is a close 
approximation of another. 
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The stability study will be carried out using Eqs. (7.61)-(7.64). In these dimensionless equations, the 
independent mesh variables at any mesh point (j, n) G fi are z/", (vs)™, (z'ss)", and (i'xxx)^- Hereafter the 
round-off errors associated with these variables are denoted by 5{vx)f , and respec- 

tively. Because the magnitude of the round-off error of a quantity with a smaller magnitude 
generally is not necessary smaller than that of a quantity with a larger magnitude, in contrast 
to Eq. (7.69), we assume that 


Wn)?| « \S(r 2 )]\ « \8v?\ « \5(u^\ « \S(u^)]\ « |<fe)?l (j,n) G fi (7.76) 

Hereafter (i) the symbol implies that the quantities on both sides of this symbol have the same order 
of magnitude, and (ii) 6(t i)" and S(t 2 )" denote the round-off errors of (n)" and ( 72 )", respectively. Also 

because round-off errors could vary erratically from one mesh point to another, in contrast to 
Eq. (7.72), in general, we have 

IK+V2 - 'K vl * K ± "i% a | « l*(n)?| « l*fo)?l (J,n) G H (7.77) 

Let x^, £ = 1,2, 3, 4, be independent mesh variables at a mesh point. Let Sx£, i = 1,2, 3, 4 be their 
round-off errors, respectively. Then because each 5xe is a very small perturbation of xe, the round-off errors 
of a differentiable function fix i,X 2 ,* 3 ,* 4 ) can be approximated by 



(7.78) 


As an example, let 

def n— 1/2 def , \n— 1/2 def , \n— 1/2 i def , \n— 1/2 / • \ ^- r\ (n nc\\ 

Xi = v j+ i /2 , *2 = {Vx) j+l / 2 i *3 = (%i) j+1/2 i and *4 = {Vxxx) j+1/2 \Ji n ) G S2 (7.79) 

Then the expression on the right side of Eq. (7.61) associated with the mesh point (j + 1/2, n — 1/2) can be 
expressed as g{x i,X 2 ,* 3 ,* 4 ) where 


g{x I,x 2 ,x 3 ,x 4 ) = f i(l - y)xi - i[l - (xi) 2 ]x 2 - ^ [3(xi) 2 - l] (1 - 2x 2 )(x 2 ) 2 


+ y [! - (*i) 2 ] (1 - 6x 2 )x 3 + ^ [! - (*i) 2 ] 2 *4 


Eq. (7.80) implies that 
dg 1 - £i 


9xi 2 


+ X 1 X 2 - 2xi(l - 2 x 2 )(x 2 ) 2 + 


[l - 3(xi) 


2x1X4 


3 (1 — 6x 2 )x 3 — ' 3 [1 - (xi) 2 


dfJ - ^ - \ [3(xr) 2 - 1] x 2 (l - 3 x 2 ) - 2xix 3 [1 - (xi) 2 ] 


9x 


2 3 

dg = xi [1 - (xi) 2 ] 
9x3 3 


- 2xiX 2 [l - (xi) 2 


and 


Note that Eq. (7.69) implies that 


dg = [1 - (*i) 2 ]~ 

9x4 6 

1*41 < 1*3 1 <| *2 1 <| *l| < 1 


(7.80) 

(7.81) 

(7.82) 

(7.83) 

(7.84) 

(7.85) 
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As such, by excluding the case 


|1 - (zi) 2 | < 1 (7.86) 

one concludes that the expression on the right side of each of Eqs. (7.81)-(7.83) is dominated by the first 
term which involves only the variable x\. Moreover, because Eq. (7.76) implies that 

\5x4,\ « |<fa 3 | \5x 2 \ « |foi| (7.87) 

With the aid of Eqs. (7.78) and (7.79), one concludes that 


6 g(xi,x 2 ,x 3 ,x 4 ) = ^ ai )&ri 




i - (*.) a i + + lh^id_ fa4 

3 6 

, ( i -- 2 ) 2 ,.. r -1/2 


2 

1 -v 2 


8v xx + 8v xxx 

3 6 


3+ 1/2 


Next, by combining Eqs. (7.75) and (7.88), one arrives at the conclusion 


Sg(xi, X 2 , X 3 , X 4 ) = s ^j+i /2 


P(1 - v 2 ) 


0 

21 


( L - V \ n 1/2 

t — 2 — A- ^W /2 


J 0 


<r/ \n— 1/2 

0\ u xx) j+1/2 


(1-D 2 ) 


2 \ 2 ^ 


6 


J i 


Cf ^n— 1/2 

0 \yxxx) j + 1/2 


Note that, according to Eqs. (7.75) and (7.79), Eq. (7.86) <t=> 

(*;) 2 = i 


(7.8 


(7.89) 


(7.90) 


As such the validity of Eq. (7.89) has been established except the case Eq. (7.90). 

By excluding the case Eq. (7.90) and using similar arguments, one can also derive the counterpart of 
Eq. (7.89) for the expression on the right side of Eq. (7.61) associated with the mesh point ( j — 1/2, n — 1/2). 
After combining Eq. (7.89) and its counterpart, one concludes that, excluding the case Eq. (7.90), 




. /I — V\ n n- 1/2 . /l + Z'V" n- 1/2 


1/2 £( \n— 1/2 

V 2 )j n S ^+V 2 ^(^)/-l/2 


1 - 1/2 


£>(1 - i> 2 )l ra ‘ 

3 J,- L 


1 - 1/2 


-1/2' 

1 

r(l-z> 2 ) 2 i 

n ■ 

!/2_ 

1 

L 6 J 

3 - 


(7.91) 


C/ \n— 1/2 _ r/ \n— 1/2 

0 \yxxx )j+\/2 0\yxxx) j_\j 2 


The proof of Eq. (7.91) is completed if one can show that it is valid even for the case Eq. (7.90). Note 
that Eq. (7.90) <t=> either 

(7.92) 

(7.93) 


*" = 1 


or 


= -! 


but not both. On the other hand, it can be shown easily that (i) 


for the case Eq. (7.92), and (ii) 


= * T /// 




(7.94) 

(7.95) 
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for the case Eq. (7.93). By comparing Eq. (7.91) with Eqs. (7.94) and (7.95), it is clear that Eq. (7.91) is 
also valid for the case Eq. (7.90). QED. 

Except for a possible complication to be discussed later, in a similar manner one can derive from 
Eqs. (7.62)-(7.64) the following results: 


6(^)1 = 


1 


L2(1 + t 2 )J, 

r 1 - t 2 - 2z> 


8v 


o+ 1/2 


8v r ‘. 


-l/2\ 

T — 75 + 2 ir 

1/2 ) 

V 2(1 + 75 ) J 


y- 1/2 

. S ( v s)j +1 /2 


x , in- 1/2 
°\ V x) j - 1/2 


(1 - r 2 + 2 v) 2 (1 + r 2 ) 


2 


L 2(1 + r 2 ) J j 
(1 - r 2 - 2z>) 2 (1 + r 2 ) 2 ■■ 


4(1+75) 12(1 + Tl) J j 

in-l/2 


r/ in— 1/2 

o( l/ a:a:) :; - + i/ 2 


.5(^)”r 1 1 /2 


4(1 + 75) 12(1 + n)Jj 

(1 + t 2 ) 2 (1 — ri + 2£) (1 — r 2 + 2j/) 3 -' 


(7.96) 


12(1 + 71) 12(1 + r 2 ) 

(1 + t 2 ) 2 (1 — 7 i — 2 z >) (1 — r 2 — 2£>) 3 - 1 




1 


12(1 + 7l) 

t-1/2 


12(1 + 75) 
1 - 1/2 


V/ \n— 1/2 

‘c/ ,n— 1/2 

e [yxxx ) j—i/2 


sm] = %\ (i - + (i + a - - 6(v sss )]:% 

and 


1 - 1 / 2 ' 


(7.97) 


fe*)" = 


r 1 1 

n ■ 

-2(1 + Tl) - 

3 - 


<■[ in— 1/2 c / in— 1 / j 

0(^x) j+ i/ 2 - <H I/ Ss) J --l/ 2 


n- 1/2 


1 — Ti + 2f> 


L 2(1 + Ti) 


r/ in— 1/2 


1 — 71 — 2z> 


(7.98) 


L 2(1 + Tl) 


r/ in— 1/2 

C l+xxx ) j — 1/2 


In deriving Eqs. (7.96) and (7.98) from Eqs. (7.62) and (7.64), the reader should bear in mind the 
simplifying convention established earlier for the parameters T\ and r 2 . As an example, we have 


n— 1/2 „ n-1 / 2 

def j±l/2 


[2(1 + r 2 ) J j±i /2 2 [l + ( 75 )"] 

Thus, by using Eqs. (7.75) and (7.78), Eq. (7.99) implies that 


(7.99) 


V 

n- 1/2 

1 

[ 2(1 + r 2 ) J 

1 + 1/2 

[ 2(1 + r 2 ) J 


5v n ; 


- 1/2 


2(1 + 12) 2 J 


S(r 2 )) 


(7.100) 


In turn, Eq. (7.100) implies that 

5 


V 

n-l/2 

-8 

1+1/2 

V 

n- 1/2 

1 

n 

( 

[ 2(1 + r 2 ) J 

[ 2(1 + r 2 ) J 

J — 1/2 

[ 2(1 + r 2 ) J 

. V 

1 


r j- n— 1/2 r n-1/2^ 

^+1/2 - ^-1/2 . 


(7.101) 


The expression on the right side of Eq. (7.101) is the first term on the right side of Eq. (7.96). 

Eqs. (7.91) and (7.96)-(7.98) can be cast into the matrix form 

Sq(j,n)= [Q-(ti, r 2 , £■)] " 8q(j + 1/2, n — 1/2) + [Q+(ri, t 2 , j/)] " 8q(j — 1/2, n — 1/2) (j,n)<E O (7.102) 
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where 


and 


with 


Sq(j,n) d A f 


( K \ 

5(l'xx) r j 


[Q±(ti,t 2 , J>)]" = f 


def 

9 n = 


1 ± V „± def , ( 

— > 912 = ±- 


(j, n) € Cl 


\6(u 

XXX ) 

U 

3 7 


/9u 

9?2 

9i*3 

9*4 \ 

9*1 

922 

923 

9m 

0 

0 

9*3 

9m 

l 0 

0 
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qtiJ 



v+ def 

ha = 

T— 


(j, n) € Cl 


„± def (1 - £ 2 ) 2 

914 = T t 


(7.103) 


(7.104) 


(7.105) 


~± def Tl def (1 ^2 “F 27') def .(It 72) 3 (1 7"2 “F 27)^ 

921 “ 2(1 + t 2 )’ 922 = 2(1 + r 2 ) ’ 923 = 12(1 + n) T 4(1 + r 2 ) 

,± def (1 + ^~2 ) 2 ( 1 ~ Tl T 27) _ (1 - T 2 =F 27) 3 
924 ~~ 12(1 + Tl) 12(1 + r 2 ) 


def 1 i ^ def , (1 7 ^) def “Fl def (1 'Tl “F 27 ) 

933 = "" 2“ ’ 934 “ 2 ’ 943 = 2(1 + n) ’ 944 ~ 2(1 + n) 


(7.106) 


(7.107) 


A comparison of Eqs. (7.102)-(7.107) with Eqs. (6.7), (5.54), and (6.8)-(6.11) reveals a striking similarity 
in the structures of these two sets of equations. On the other hand, it was explained in Sec. 5.4 that: (i) 
computational instability would occur if round-off errors are amplified without bound as they propagate down 
the subsequent time levels; and (ii) for linear schemes like the a-e-4 and c-t -4 schemes, the time evolution of 
the mesh variables and that of the associated round-off errors are governed by the same system of equations. 
As a result of the facts stated above, one now arrives at the conclusion that the stability conditions of 
the nonlinear inviscid Burgers c-t -4 scheme can be obtained from those of the linear c-t scheme 
by replacing the constant coefficients ;/, ti, and t 2 with the mesh-point dependent coefficients 
7 ”, (ti)!?, and (t 2 )™, respectively. 

7.3. The inviscid Burgers c-t *-4 schemes 

Let fti(s) and ft. 2 (s) be any two functions that meet the same requirements imposed on the function 
h(s) introduced in Sec. 3 . 2 . Then, because of the close relations among the c-t, c-t -4 and current schemes 
established above, here an “ideal” inviscid Burgers c-t *-4 scheme is defined by an extension of Eq. ( 6 . 21 ), 
i.e., 

(n)" = /ir((j>") 2 ) and (t 2 )” = h 2 {(i>?) 2 ) (l*"l < 1 ) ( 7 - 108 ) 

According to Eq. (3.22), one can assume hi(s) = 7 2 (s) = yfs. As such a special inviscid Burgers c-t *-4 
scheme can be defined by the relation 


(n)? = ( r 2 )] = |^| (1^1 < l) 


(7.109) 
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For robust real-world applications, we may also introduce a special inviscid Burgers c-r-4 scheme which 
is defined by an extension of Eq. (6.23), i.e., 


(n)" = Pi \i>j\ and fo)? = 02 W\ (Pi > l; p 2 > l; \i>?\ < l) (7.110) 

where Pi > 1 and P 2 > 1 are two adjustable parameters. In the rest of this section, we consider only the 
special schemes which satisfy at least one of Eqs. (7.59), (7.60), and (7.108)-(7.110). As a result of Eq. (7.43), 
for these schemes, 

( ti )" — > 0 and (r 2 )" — » 0 as At — > 0 (7.111) 

7.4. Shock capturing schemes 

As explained in Sec. 2.3, the current Burgers solver needed to be modified so that it is capable of 
resolving shocks. Specifically, we need to turn the expression on the right side of Eq. (7.52) and/or that of 
Eq. (7.58) into proper weighted averages discussed in Sec. 2.3 and Sec. 4. 

First we consider Eq. (7.52). Let 


,(p + ) = (f)‘ 


r(P + ) and u S i(P ) d = 


def ^A Xy 


PP~) 


(7.112) 


where u xx (P + ) and u xx (P~) are defined by Eqs. (7.47), (7.48), (7.50), and (7.51). Then it can be shown 
that 

u xx (P + ) = [(1 - 6 v x )u xx - (1 - n + 2 v)u xxx \ ” /+ ^ (7.113) 

and 

u xx (P~) = [(1 - 6 v x )u xx + (1 - ri - 2v)u xxx ] n j _l / / l (7.114) 

Thus Eq. (7.52) implies that 

( u xxx)j — 2 1 + (r )" ~ 2 ( Uxxx ~ P u xxx+)j (7.115) 


where 


V Uj xxx— ) j — 


def (Uxx)^ u xx(P ) 


1 + (nYi 


and ( u xxx+ ) r j' = 


def u xx(P^~) ( u xxYn 


l + (n)? 


(7.116) 


with the understanding that ( u xx )” is that defined in Eq. (7.44). Eqs. (7.115) and (7.116) can be considered 
as the substitution forms of Eqs. (3.10) and (4.1)-(4.3). 

According to Eq. (7.115), (u X xx)(j evaluated using Eq. (7.52) is the simple average of (' u xxx _)” and 
(u xxx +) r j- Moreover, by using Eqs. (7.111), (7.113)-(7.116) and 


v’j -»■ 0, (; u x ) 1 -»■ 0, , (v xx Y- -> 0, and (v xxx ) T t 


as At — > 0 


(7.117) 


(which follows from Eq. (7.43)), one can show that 


lim (u xxx —)j — lim (Uxxx+Y — lim (u xxx Y 
At— »0 J Af->0 J Af— »0 J 



XX ^ 


XXX 


xn-1/2 

'j+ 1/2 


(u xx + u 


XXX 


\n- 1/2 


(7.118) 


Because (■ u xxx -)" and (u xxx +)(j satisfy Eq. (7.118), the simple average in Eq. (7.115) can be turned into 
an weighted average using any of the techniques described in Sec. 4. As an example, consider the technique 
presented in Sec. 4.4. The current counterparts of Eqs. (4.22)-(4.24) are 


[M 


1 n def 
± ; = 


I (^xxxzt ) 7 


n I ai 

ij 


min{ | (u xxx -) r l | ai , | (u xxx+ ) r l | Ql } 


- 1 


(ai > 0 and \(u xxx -)j \ + \(u xxx+ )]\ > 0) (7.119) 
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and 


v / ~ \ in def 

,(wi)+L = 


[<».)-]," = { 


1 + 0-1 ( vi )- 


2 + oi[(r 7 i)+ + (771 )-}\ j 

i+qi(t?i)+ y 

2 + o-i[(r?i)+ + (r?i)_] j , 


iy? * 0 ) 

(7.120) 

(*? y 4 0) 

(7.121) 


respectively. Here (i) aq > 0 is a preset parameter in the order of 1, and (ii) 


{ \n 

( 01)7 — 


def (oi) 0 


((cr 1)0 > 0 and v 7 - ± 0) 


(7.122) 


with (a i) 0 being a preset parameter in the order of 1. With the above definitions, the current counterpart 
of Eq. (4.30) can be expressed as 


(«sss)j = [(tci + [(* 1 )+] '-(Uxxx+)j (7.123) 

Next consider Eq. (7.58). Note that, with the aid of Eq. (7.52) and (7.57), it can be shown that Eq. (7.58) 




_ . 4 -(C)-.*(r) _ [i + fog] , 

— oh , I- \ n.l a \Uxxx)j 


2 [l + ( T 2)j 


6 


In the following, with the aid of Eq. (7.124), (u$±)™ will be defined so that 


(' u x)j = ~(Ux-+U s+ )] 


and 


lim (Un-)] = lim (u $+ )” = lim (u s )j = - 
J Ai— >0 J Ai— >0 J 2 


Ai— >0 


Ai— >0 


n- 1/2 


U - U S + -Uxx I -(» + % + r«xx 

d '7+1/2 V d + 7 — 1/2 


n- 1/2' 


(7.124) 


(7.125) 


(7.126) 


Note that the validity of the last equality sign follows directly from Eqs. (7.58), (7.111), and (7.117). 

To facilitate the development, note that here (i) P + and P~ are the points depicted in Fig. 3, except 
that the parameter r shown in the figure be replaced by the mesh-point dependent parameter ( 72 )"; and 
(ii) u‘ l *(P + ) and u A *{P~) defined in Eq. (7.56) are 4th-orcler Taylor’s approximations of the values of u at 
points P + and P ~ , respectively. Alternatively, these two values can also be approximated by 


u 4 *(p±) a* u * (x(P ± ),t n ; j,n) 

respectively. With the aid of Eqs. (7.35), (7.42), (7.53), and (7.54), Eq. (7.127) implies that 

,. + lu r-n,v?' I 4- ^ + 1—1+ 

3 


U 


4 *(P±) = U] ± [1 + (+ 2 )"] («x)? + [l + ( 0 T2);i] («xx)? ± [l + ( 7 )j?] («xxx)? 


In turn, Eq. (7.128) implies that 


~4* 


(p+) - 


1 + M 


6 


-« 4 *(P-) , 1 + (r 2 )j 


J -L( u --) n = 1 1 _j v ( u T 

1 +(r 2 )7 2 1 1 + ( r 2 )" 0 1 ' 7 


(7.127) 


(7.128) 


(7.129) 


2 v 1 + (r 2 )” 2 

Moreover, by using Eqs. (7.44), (7.46), (7.52), (7.57), (7.58), (7.111), (7.117) and (7.128), one concludes that 


lim ft 4 *(P ± ) = lim u 4 *(P ± ) = [« + «! + + - Uxxx 

At—>0 K ’ At—>0 K ’ V 2 xx T g xxx 


x— 1/2 


j± 1/2 


(7.130) 


49 



An immediate result of Eqs. (7.129) and (7.130) is 


lim (?4 + )" = lim (4-)? 

Ai^O x+ 3 At-^0 x J 


where 


and 


(4+)" = 


def u 4 *(P+)-u " 1+4)) 


1+4)) 


( u-xx ) j 


, S„def4- U ( P “) , 1 + 4) 


(4-)? = 


i + 4)," 

As a result of Eqs. (7.132) and (7.133), one has 


-447 


u r + u ) -- “ 4, ( p+ )-“ 4 *( p -) 


(7.131) 


(7.132) 


(7.133) 


(7.134) 


Combining Eqs. (7.124), (7.131) and (7.134), one concludes that the requirements Eqs. (7.125) and (7.126) 
are met by the definition 

/ \n def f / n [l + ( ^”2 ) j ] . , n . . 

\ u x±)j — ( u x±)j g \ u xxx)j (7.135) 

Because (ux+)j‘ and (ttg _) 7 meet the condition that they have the same value in the limit of At — > 0, 
the simple average on the right side of Eq. (7.125) can also be converted into an weighted average using 
any of the techniques described in Sec. 4. As an example, consider the technique presented in Sec. 4.4. The 
current counterparts of Eqs. (7.119)-(7.121) are 


[4)4? = 


144 


n \ot2 
3 1 


j min{|(u s _)7| a2 , 14474} 

def f 1 + ct 2 (?7 2 )_ 


.44 


and 


[(*)-]”"{ 


2 + o- 2 [(r? 2 ) + + 4 )_] 7J 

def / 1 + cr 2 40 + 


2 + cr 2 [(772 )_|_ + (772)-] 


-)?l + I4+)?I > 0 ) 

(7.136) 

(4 * 0 ) 

(7.137) 

(4 ^ 0 ) 

(7.138) 


respectively. Here (i) a 2 > 0 is a preset parameter in the order of 1, and (ii) 

(( 0 - 2)0 > 0 and 4 7 ^ 0 ) 


/ \n def (0-2)0 
( 02)7 — l--.nl 


(7.139) 


with (o - 2 ) 0 being a preset parameter in the order of 1. With the above definitions, the current counterpart 
of Eq. (7.123) can be expressed as 


4)7 = [4 2 )-]”4-)7 + [4 2 )+]>47 (7.i40) 

Note that a result of Eqs. (7.125) and (7.135) is 

4)? = \ (4+ + 4-); - [1 + ( 6 T2)j?] 44? (7-i4i) 

Thus an alternative shock-capturing form of Eq. (7.125) can be obtained by replacing the simple-average 
term on the right side of Eq. (7.141) with an weighted average of (4+)7 an d (4-)?- 
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8. A sketch of multidimensional extensions 


Consider the 2D advection equation 


du 

~dt 


du du 

d X TT 1" dy TT" = 0 

dx dy 


( 8 . 1 ) 


where a x , and a y are constants. Let x% = x, X 2 = y, and X 3 = t be the coordinates of a three-dimensional 
Euclidean space E 3 . By using Gauss’ divergence theorem in the space-time £ 3 , it can be shown that Eq. (8.1) 
is the differential form of the integral conservation law 


<£ h-ds = 0 (8.2) 

JS(V) 

Here (i) S(V) is the boundary of an arbitrary space-time region V in E 3 , (ii) 


7 * def 

a = 


(a x u, a y u, u) 


(8.3) 


is a current density vector in £ 3 , and (iii) ds = dan with da and n, respectively, being the area and the 
outward unit normal of a surface element on S(V). 

In addition, by using an argument similar to that was used to obtain the conservation condition Eq. (5.9), 
it can be shown that, for a smooth solution u(x,y,t) of Eq. (8.1), we have 


where 



• ds = 0, 



• ds = 0, 


and 


® h yy ■ ds = 0 
JS(V) 


h ^ef &h ( (Pu cPu, d 2 u \ h ^ef d 2 h 

xx dx 2 \ 1 9a: 2 ’ v dx 2 ' dx 2 ) ’ xy dxdy 


and 


r def d' 2 h 

hyv ~ Q^2 


f d 2 u d 2 u d 2 u\ 


( d 2 u d 2 u d 2 u \ 
\ x dxdy ’ Uy dxdy ’ dxdy ) 


(8.4) 


(8.5) 


In a 2D CESE scheme [13], a spatial domain can be divided into triangles of arbitrary shape with the 
understanding that any two neighboring triangles share a common side. Each triangle and the associated 
spatial solution point is assigned with an integer cell number j. The x- and y- coordinates of a spatial solution 
point j are denoted by Xj and yj , respectively. In addition, a space-time solution point (j, n) is defined to 

def 

be the point with x = xj, y = yj, and t = t n = nAt. As in the ID case, each (j. n) is associated with a 
solution element SE (j,n). However, corresponding to the three sides of a triangle, each (j, n) is associated 
with three basic space-time conservation elements CE r (j, n), r = 1,2,3. Again these BCEs are constructed 
so that the boundary of each BCE can be divided into component parts with each of them belonging to a 
unique SE. The details are described in [13]. 

To construct a 4th-order scheme, for any ( x,y,t ) £ SE the dependent variable u is approximated 
by a numerical analogue of the third-order Taylor’s expansion, i.e., 


u*(x,y,t;j,n) d = u] + {u x )](x - xj) + {u v )]{y - yj) + - t n ) + ^(u xx )](x - Xj ) 2 

+ \( u vv)j(y - Vi ) 2 + \{ u ttYj{t - t n ) 2 + (u X y)^(x - Xj)(y - Vj ) + (u xt )]{x - Xj)(t - t n ) 
+ ( Uyt)j(y - Vj){t - t n ) + ^{u xxx )](x - Xj ) 3 + ^(u VV y)](y - yj ) 3 + i(titit)"(i - t n ) 3 
+ ~^{ u xxy)j {% ~ x j) ( y — l Ij) + 2 ( Uxyy )j ( x ~~ x i)(y — Vj) d" 2 ( u xxt)j { x — x j) (t — t ) 

+ \{ u vvt )"(y - Vj) 2 (t - t n ) + ^(u xtt )j{x - Xj)(t - t n ) 2 + ^(u ytt )j(y - y 3 )(t - t n ) 2 
+ ( u xy t)j ( x — Xj)(y — Vj)(t — t ) 
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Let u = u*{x,y,t\ j,n) satisfy Eq. (8.1) uniformly in SE i.e. 


du*(x,y,t\j,n ) du*(x,y,t;j,n) du*(x,y,t;j,n) _ 


dt 


Then it can be shown that 


dx 


dy 


= 0 


(x, y, t) € SE (j,n) 


and 


(^t) j — (®x‘U'x T dy Uy ) j 

i^xt) j — fax'Uxx T dy'U'xy) j 

( u yt)j — — ( a a ;U X y + a y u yy) j 

(utt ) j = (.Q j x‘U'xt T a y u yt ) = [(^x) ^xx T ^^x^y'Uxy T i^y) '^-yy\ j 
(^xxt) j i^x^xxx T QyHxxy) j 

(U xy t)j = -( Q j x'U j xxy “1“ Q j y r U j xyy) j 
(' u yyt)j = ~{ a xUxyy + a y u yyy)j 

( l U , xtt)j ( Q'x'U'xxt T Uy'^xyt ) j — [(^x) ‘U’xxx T ‘^‘O’x^y ^xxy T ( fly ) U-xyy\ j 
(■ u ytt )] = ~{a x'U'xyt H - Q'y^yyt) j — [(^a?) 'U j xxy “I” < 2‘Q’xQ j y'U j xyy “1“ if^y) ^yyy\ j 

( u ttt)j = ( a x u xtt H - a y u ytt)j = [(^tc) u xxx H" 3(d x) dyllxxy T H xyy T (&y) 


l yyy J 


Substituting Eqs. (8.8)-(8.17) into Eq. (8.6), one has 


(8.7) 

( 8 . 8 ) 

(8.9) 

(8.10) 

( 8 . 11 ) 

(8.12) 

(8.13) 

(8.14) 

(8.15) 

(8.16) 

(8.17) 


u*(x,y,t;j,n) = u] + (u x )j[(x - xj) - a x {t-t n )\ + (u v )][(y- yj) - a v (t-t n )\ 

+ ^(u X x)j[{x - Xj) -a x (t-t n )] 2 + ^(u yv )][(y - yj ) - a v (t - t n )] 2 

1 3 

+ (u xy )j [{x - Xj) - a x (t - t n )\ [(y - yj) - a y (t - t n )] + -(u xxx )] [{x - Xj) - a x (t - t n )\ ( 8 . 18 ) 

^(«wj)i[(s - Vo) ~ a y(* - i ")] 3 + \( u xx V )j[(x - Xj) - a x (t - t n )Y[{y - Vj ) - a y (t - t n )] 

+ \( u xyy)j[(x -Xj)- a x {t - t n )\ [( y -Vj)- a y {t - t n )\ 2 

Thus u*(x, y, t ; j, n) is a 3rd-order polynomial function of [(a: — Xj) — a x (t — t n )] and [(y — yj) — a y (t — t n )] 
specified by the coefficients u ™, (u x )", (%)”, («xx)", {u vy )j, (uxy)", ( u XX x)j , (u yyy )j, ( u X xy)j , and (u xyy )™. 
For a reason to be given later, in the current extension, only eight of these ten coefficients are considered 
as the independent mesh variables at (j, n) while two of them (i.e., ( u xxy )" and (u-xyy)") are considered as 
functions of another four independent mesh variables at ( j,n ). In other words, there are twelve independent 
mesh variables at each solution point ( j,n ). 

To construct the o-4 scheme for Eq. (8.1), we will impose the following three conservation conditions: 


(f h*-ds = 0, r = 1, 2, 3 (8.19) 

JS(CE r {j,n )) 

where 

h*(x, y, t ; j, n) = [a x u*(x,y,t-,j,n),a y u*(x,y,t;j,n),u*(x,y,t-,j,nj\ (8.20) 

Eqs. (8.19) and (8.20) are the numerical analogues of Eqs. (8.2) and (8.3), respectively. 
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Note that Eq. (8.18) implies that 


d 2 u*(x,y,t\j,n ) 
dx 2 

d 2 ii*(x, y, t ; j, n) 
dxdy 


and 


d 2 u*(x,y,t\j,n) 
dy 2 


= {u xx )j + (u xxx )j [(a: - Xj) - a x (t - t n )] + ( u xxy )j [(y - y 3 ) - a y (t - t n )\ (8.21) 


= ( u xy )" + ( u xxy )" [(a; - Sj) - a x {t - t n )\ + (u xyv )™ [(y - yj) - a y (t - t n )] (8.22) 


— ( U yy)j + ( u xyy)j ( x x j ) a x(t t ) + ( Uyyy)j [(V Vj) 0, y (t t ) (8.23) 


Thus, at each (j, n), imposition of conservation conditions Eq. (8.19) and the numerical analogue of Eq. (8.4), 
i.e., 


(f h* ■ ds = 0, 

Js(CE r (j,n)) 

r* = 1,2,3 

(8.24) 

f h* y -ds = 0, 

Js{CE r (j,n)) 

r = 1,2,3 

(8.25) 

K y -ds = 0 , 

JS(CE r (J,n)) 

r — 1,2,3 

(8.26) 


and 


would result in twelve conditions for the ten coefficients that appear on the right side of Eq. (8.18) if one 
assumes that (i) 


and 


and (ii) 


and 


K x {x,y,t;j,n) d = [a x u* xx (x,y,t\j,n), a v u* xx {x,y,t\j,n), u* xx (x,y,t;j,n)\ 
h xy {x, y, t; j, n) = f [a x u* xy {x,y,t-,j,n), a v u* xy {x,y,t-,j,n), u* xy (x,y,t;j,n )] 

K y (x,y,t-,j,n) = f [a x u* yy (x,y,t\j,n), a y u* yy (x,y,t-,j,n), u* yy (x,y,t;j,n)\ 


u* xx {x,y,t\j,n) 
u xy (x,y,t;j,n ) 


d 2 u* (x,y,t\ j, n) 
dx 2 

d 2 u*(x,y,t;j,n) 

dxdy 


u y V {x,y,t;j, n) 


d 2 u* (x,y,t\ j, n) 
dy 2 


To overcome the above problem of over determinacy, we replace Eqs. (8.30)-(8.32) with 


Kx(x, y, t\ j , n) = ( u xx )!? + {u xxx ) n j [(a; - xj) - a x (t - t n )] + (■ u xxy )? [(y - yj) - a y (t - t n )] 

U xy{Xl U: jl n ) = ( U xy)j + (' U xyx )j [(a; — Xj) — <X x (t — t )] + (u a ;yy ) j \_(y — yj) ~ a y{t — t )] 

u yy ( x , y , t\ j , n) = (u y y)j + (• u yyx )j [(a; — Xj) — a x (t — t )] + (■ u yyy )j [(y — yj) — a y (t — t )] 


(8.27) 

(8.28) 

(8.29) 

(8.30) 

(8.31) 

(8.32) 

(8.33) 

(8.34) 

(8.35) 
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respectively. Note that: (i) the expression on the right sides of Eq. (8.33) is identical to that of Eq. (8.21) 
except that the symbol (w xxy )" in Eq. (8.21) is replaced by ( u xxy )" in Eq. (8.33); (ii) the expression on 
the right side of Eq. (8.34) is identical to that of Eq. (8.22) except that the symbols ( u xxy )" and ( u xyy )" 
in Eq. (8.22) are replaced by ( u xyx )" and (« xyy )", respectively; and (iii) the expression on the right side 
of Eq. (8.35) is identical to that of Eq. (8.23) except that the symbol (u xyy )" in Eq. (8.23) is replaced by 
( u yyx )”. In spite of the facts that (i) (u xxy )”, (tt xxy )", and ( u xyx )" are the “numerical analogues” of the 
value of the same mixed spatial derivative d 3 u/dx 2 dy at (j, n); and (ii) (u xyy )™, ( u xyy )”, and ( u yyx )” are the 
“numerical analogues” of the value of the same mixed spatial derivative d 3 u/dxdy 2 at (j, n), here (i) (u xxy )", 
(i( XX9 )“, and ( u xyx )j are treated as different numerical parameters and thus can assume different 
values; and (ii) (u xyy )”, ( u xyy )"* and ( u yyx )" are also treated as different numerical parameters 
and thus can assume different values. In fact, (u XXJ/ )™, (« xyx )", ( w xyy )", and (u yyx )” are considered 
as independent mesh variables at (j, n) while (i) (u xxy )" is considered as a function of (u xxy )™ 
and (u xyx )"; and (ii) (it xyy )" is considered as a function of (w xyy )" and (u yyx )™. In fact, one may 
assume that 

[u xxy ) j = d U’xyx) j and {u xyy ) j = — (u xyy + u yy x)j (8.36) 


As such, at each (j, n), there are twelve independent mesh variables, i.e. , u”, (u x )", (« y )", (u xx )", 
i'^xxx} j , {v j xxy) j ■ i.'U'xy) j ; i^xyx)j , (^xyy)j , ( ^yy)j 1 (^yyx) j ■ and (u yyy jy • They can be determined in terms 
of the known mesh variables by imposing the twelve conservations conditions given in Eqs. (8.19) and (8.24)- 
(8.26). In fact, with the aid of Eqs. (8.27) and (8.33), the three mesh variables (u xx )", (u xxx )", and (u xxy )" 
can be determined in terms of the known mesh variables by imposing the three conditions given in Eq. (8.24). 
Similarly, (u xy )", (u xj/x )”, (u xyy )™ can be solved by using Eq. (8.25) while (u yy )j, (u yyx )™, (u yyy )j can be 
solved by using Eq. (8.26). After the nine 2nd-order and 3rd-order independent mesh variables at (_), n) 
are determined, the last three independent mesh variables, i.e., u", (u x )", and (w y )" can be solved using 
Eq. (8.19). Thus the a-4 scheme for Eq. (8.1) can be constructed easily. Moreover, with the aid of Eqs. (8.33)- 
(8.36), the dissipative extensions of the 2D a-4 scheme can also be constructed easily using well-established 
CESE techniques. 

We conclude this section with the following remarks: 

(a) The identities 


9 2 _ d 2 fdu\ 

dx 2 \ <9y y \ dx J 


d 2 rdu\ 

dydx \ dx J 


(8.37) 


and 


d 2 / du\ 
dy 2 \ dx J 


d 2 f du\ 
dxdy \dy J 


d 2 f du\ 
dydx \dy ) 


(8.38) 


are derived from a limiting process in which ax, Ay — > 0. Because Ax and Ay are finite in a numerical 
simulation, logically the numerical analogues of the third-order mixed spatial derivatives are not required 
to satisfy the numerical analogues of Eqs. (8.37) and (8.38). 

(b) Because ( u xxy )" and (u XJ/y )" appear as 3rd-order terms on the right side of Eq. (8.18), the 4th-order 
accuracy requirement of u*{x,t\ j,n) demands only that (u xxy )” and (u xyy )" be lst-order in accuracy, 
i.e., the same orders of accuracy achieved by (w x )” and (u y )" in a 2nd-order 2D CESE scheme. This 
rather low requirement for accuracy further justifies the above introduction of (u xxy )", (u xyx )”, (u xyy )", 
and (u yyx )j. 

(c) For the 3D CESE method, tetrahedrons are the most natural building blocks for structured or unstruc- 
tured spatial meshes. As it turns out, the high order CESE framework described here can be naturally 
extended in a 3D space using spatial meshes built from tetrahedrons. 


9. Numerical results 
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To assess the stability and accuracy of the a-e-4 and c-r-4 schemes, consider the model problem defined 
by the PDE Eq. (2.1) and the exact solution 


u = u e (x, t) = sin [27r(x — t )] 


(9.1) 


We have 


a = X=T = 1 


where A = wavelength and T = period. Let (i) 

, . def du e (x,t) det d 2 U e (x,t) , def d 3 U e (x,t) 

U xe {X,t) = — , U xxe {x,t) = ^2 1 and U xxxe (x,t) = ^3 

and (ii) the spatial domain of unit length be divided into K uniform intervals. Thus 

Ax = 1 /K, At = pax and t = nAt 


(9.2) 


(9.3) 


(9.4) 


where n = number of marching cycles (each cycle is formed by two consecutive marching steps), and t = 
total marching time. 

In Tables 1-14, the numerical errors of several computations using the a-e-4 and c-r-4 schemes are 
presented in terms of the L 2 error norms 


E(K, n, v) 


N 


1 

K 


K 


5>" 


U e (Xj,t n )\ 


2 


(9.5) 


and 


E X (K, n, p) d = 


\ 


E XX (K , n, v) d = 


N 


E XXX {K , n, v) = f 


\ 


JT ±[(u x )? - U xe {Xj,t n )] 2 
3 = i 


1 x V 

y^t^xx)? - u xxe (xj,t n )] 2 
A i=i 




- 


=(*i »*")] 2 


(9.6) 

(9.7) 

(9.8) 


Note that {u x )™, (u x a,)", and (it Jxa; )" which appear in Eqs. (9.5)-(9.8) are not the normalized mesh variables 
defined in Eq. (7.42). 

The numerical errors of the a-e-4 and c-r-4 schemes with periodic boundary conditions are shown in 
Tables 1-7. In all simulations, the values of n are adjusted so that t = 5.4 for all the results presented. 

The results with v = 1 generated using the a-e-4 scheme are given in Table 1 (ei = ti = 0) and Table 2 
(ei = 62 = 0.5). For these cases where the value of v is right at the stability boundary, the errors for «" and 
(■ u xx )" are machine zero. On the other hand, the values of E x and E xxx are reduced by factors of about 16 
and 4, respectively, as both I\ and n double their values. Thus, for both the cases (a) v = 1 and ei = €2 = 0; 
and (b) v = 1 and e\ = €2 = 0.5, the a-e-4 scheme is 4th-orcler and 2nd-order in accuracy for (u x )" and 
(u xxx )™, respectively. 

The results with v — ei = ti = 0.5 generated using the a-e-4 scheme are shown in Table 3. It is seen 
that, for this case, the a-e-4 scheme is 4th-order in accuracy for both u" and (mj)" while it is 2nd-order in 
accuracy for both («„)" and (■ u xxx )^ . 

The results with v = 0.05 and ei = €2 = 0.5 generated using the a-e-4 scheme are shown in Table 4. 
It is seen that, for this case with a very small value of p, the a-e-4 scheme has an accuracy higher than 
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4th-order for both it” and (it x )” when the meshes used are relatively coarse. For the coarse mesh cases, it 
also has an accuracy higher than 2nd-order for both ( u xx )” and (u xxx )”. According to the discussions given 
in Sec. 2 and 5, the special scheme with e\ = £2 = 0.5 (referred to as the c scheme) has a huge advantage 
in computational efficiency over other a-e-4 schemes. As such the c scheme is emphasized in our numerical 
study of the a-e-4 scheme. 

The numerical results of the three cases with (a) v = n = t 2 = 1, (b) v = n = 72 = 0.5, and (c) 
v = n = t 2 = 0.05, are presented in Tables 5-7, respectively. According to a discussion given in Sec. 6, 
all these cases are the so-called c-r*-4 schemes which form a special class of the c-r-4 schemes and are 
computational efficient and also Courant number insensitive. For case (a), the errors for a" and (u xx )” are 
again machine zero while it is 4th-order and 2nd-orcler in accuracy for (uj)" and (u xxx )”, respectively. For 
case (b), again the scheme is 4th-order in accuracy for a". Also it is 3rd-order, 2nd-order and lst-order in 
accuracy for (u x )”, (tt xx )", and (it xxx )” when the meshes used become finer and finer. In fact, for this case 
we have (i) E = 2.45 x 10 -8 , E x = 2.14 x 10 -6 , E xx = 7.31 x 10 -3 , and E xxx = 0.299 when K = 200 and 
n = 2160; and (ii) E = 1.53 x 10 -9 , E x = 2.76 x 10 -7 , E xx = 1.83 x 10 -3 , and E xxx = 0.161 when K = 400 
and n = 4320. For case (c), as the meshes used become finer and finer, the scheme converges to 4th-order, 
3rcl-order, 2nd-order, and lst-order in accuracy for it", (it x )”, (it xx )", and (u xxx )", respectively. 

According to discussions presented in Sec. 2, 3, 5, and 6, in the case where \u\ <C 1, the c-r*-4 schemes 
should be much more accurate than the 4th-order c scheme. This conclusion is indeed consistent with the 
numerical results presented in Tables 4 and 7. 

The results presented in Tables 8-14 are for the cases which are identical to those associated with the 
results presented in Tables 1- 7 except that the periodic boundary conditions used in the latter cases are 
replaced in the former cases by the boundary conditions in which (i) the exact solution Eq. (9.1) is specified 
at the left spatial boundary x = 0, and (ii) the 4th-order non-reflecting boundary conditions Eqs. (5.83), 
(5.86), (5.88), and (5.89) are imposed on the right spatial boundary x = 1. As shown by the results presented 
in Tables 8-14, this changes of boundary conditions has no impact on the orders of accuracy for the a-e-4 and 
c-t*-4 schemes. In fact, this is true even for the numerical results associated with the right spatial boundary 
x = 1. 

The last numerical results to be presented is a shock solution for the problem defined by Eq. (7.1) and 
the initial condition 

( 1 if x < 0 

u = l (t = 0) (9.9) 

l 0 if x > 0 

The exact solution of this problem is 


u e (x,t) 


1 if x < t/2 
0 if x > t/2 


(9.10) 


The numerical solution is generated using a scheme described in Sec. 7. Specifically, we assume 
Eqs. (7.44), (7.46), (7.110), (7.119)-(7.123) and the alternative shock-capturing form of Eq. (7.125) described 
following Eqs. (7.141) (the same procedure Eqs. (7.136)-(7.140) is used to obtain the weighted average of 
( u 'x+)'j an< i ( u $-)j )• We also assume that 

/3i = #2 = = « 2 = 2 and (cti) 0 = ( a 2 )o = 1 (9.11) 


The spatial computational domain is defined by —0.505 < x < 0.505 and divided into 101 uniform 
interval, i.e ax = 0.01. We also assume that At = 0.009 and n = 50, i.e., t = nAt = 0.45. According to 
Eq. (9.10), 


u e (x, 0.45) 


1 if £ < 0.225 
0 if £ > 0.225 


(9.12) 
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The numerical solution is shown against the exact solution in Fig. 4. Note that even though the maximal 
CFL number = 0.9 for the exact solution, numerically, the maximal CFL number reached during the 
simulation is 0.967. 

10. Conclusions and discussions 

A thorough and rigorous discussion of a new approach for constructing highly stable high-order CESE 
schemes has been presented. Each of the linear and nonlinear 4th-order schemes described in Sec. 5-7 has 
the same stencil and same stability conditions of its 2nd-order version and still retains all other advantages of 
the latter scheme. In particular, the c-r-4 scheme and its extensions for solving the inviscid Burgers equation 
Eq. (7.1) are completely explicit, i.e., implementation of these schemes does not require inverting systems 
of linear /nonlinear algebraic equations locally or globally. As such the new approach can overcome several 
common shortcomings of traditional high order schemes that were listed in the abstract. 

Because the new 4th-order schemes are built from their 2nd-order versions, in addition to providing a 
review of the CESE core ideas and advantages in Sec. 1, a rather extensive review of several key 2nd-order 
CESE schemes were provided in Sec. 2-4. The special CESE techniques (e.g., the advanced weighted-average 
techniques described in Sec. 4) and the stability conditions presented in this review form a foundation for 
the 4th-order schemes discussed in Sec. 5-7. 

To explain clearly how a 4th-order scheme can have the same stencil and same stability conditions of 
the 2nd-order scheme from which it is built, the conceptual leap that underlines the current development 
is clearly explained early in Sec. 5. This key idea is then thoroughly validated in the rest of the paper. 
The validation includes rigorous discussions of the stability conditions of the new 4th-order schemes. In 
particular, in Sec. 7.2, we present a stability study of a nonlinear scheme based on a rather exhaustive 
analysis of the propagation of round-off errors. As a result of this study, the reader can understand why the 
underlying idea is applicable even in a nonlinear case. 

To complement the new 4th-order schemes, a set of simple 4th-order non-reflecting boundary conditions 
for linear schemes is presented in Sec. 5.5. Its potency is validated in Sec. 9. The extensions of these 
boundary conditions for nonlinear schemes are under development. 

By introducing new independent mesh variables (e.g., those appear in Eqs. (8.33)-(8.35)) and the 
concept of dependent mesh variables (e.g., those defined in Eq. (8.36)), it is shown in Sec. 8 that the 
framework described in Sec. 5-7 can be extended to construct multidimensional high order schemes. 

Accuracy and stability of the new 4th-order CESE schemes are rigorously assessed using the numerical 
results presented in Sec. 9. It is shown that the theoretical predictions made in Sec. 5-7 are consistent with 
the numerical results. 

Even though they are not presented in the paper, from the details provided here, in a straightforward 
manner, the reader can build from a given 2nd-order CESE scheme its 6th-, 8th-,... order versions which 
have the same stencil and same stability conditions of the 2nd-order scheme. 

Finally note that the current approach can also be used to construct high-orcler schemes from classical 
2nd-order schemes. In fact, by using a “regular” space-time mesh in which (i) the set of mesh points used 
include all (j, n ) with j, n = 0, ±1, ±2, . . ., and (ii) Xj = j ax and t n = riAt, a 3rd-orcler hybrid scheme has 
been built from the Lax-Wendroff and the second-order a scheme. To construct this scheme, Eq. (5.4) is 
replaced by 

u*(x,t; j, n) = uj + (u x )j [(a; - Xj) -a(t- t n )] + t(u xx )" [(a; - Xj) - a(t - t n )\ 2 (10.1) 

Also, to account for the new mesh structure, we have 

V = a ^, («*)?=* ^(«*)"> and ( u$2 )« d A f (^) 2 ( Ura )«, j, n = 0 ± 1, ±2, . . . (10.2) 

Then the new hybrid scheme is formed by 

/ \n v{y T 1)/ \n— 1 | /-i \n— 1 . v(y 1) / \n— 1 /in q\ 

yUxxjj — 2 \U j xx)j — i T (1 v )\U>xx) j T 2 \U'xx)j - |-i (10.3) 
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(10.4) 


and 


where 


and 


(«*)" = \ [(s+)"+i - (Op 1 ] - t 


\n— 1 


\n-ll 


/ _ \n def 


(O? d = 


, , 2(1 + v + 

U - (1 + v)Ux H 


J 3 


x 2(1 — ^ + i/ 2 ) 

u + (1 - v)u s H 


f («**)? 

(10.5) 

j,n = 0, ±1, ±2, . . . 

(10.6) 

j,n = 0,±1,±2, ... 

(10.7) 


J 3 


As expected, the hybrid scheme is third order in accuracy and is stable if v 2 < 1. 
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Table 1. Numerical results of the a -e -4 scheme with periodic boundary conditions 
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K=50, n=270 

K=100, n=540 
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E XX 
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E XXX 

1.09 
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Table 2. Numerical results of the a -s -4 scheme with periodic boundary conditions 
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Table 3. Numerical results of the a -s -4 scheme with periodic boundary conditions 
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Table 4. Numerical results of the a -e -4 scheme with periodic boundary conditions 
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Table 5. Numerical results of the c-r-4 scheme with periodic boundary conditions 


■- = 1.0, t = 5.4 
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Table 6. Numerical results of the c-r-4 scheme with periodic boundary conditions 
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Table 7. Numerical results of the c-r-4 scheme with periodic boundary conditions 
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t 2 =0.05, t = 5.4 
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Table 8. Numerical results of the a -s -4 scheme with non-reflecting boundary conditions 



Table 9. Numerical results of the a -s -4 scheme with non-reflecting boundary conditions 



Table 10. Numerical results of the a -e -4 scheme with non-reflecting boundary conditions 









































































































Table 11. Numerical results of the a -e -4 scheme with non-reflecting boundary conditions 



Table 12. Numerical results of the c-r-4 scheme with non-reflecting boundary conditions 
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Table 13. Numerical results of the c-r-4 scheme with non-reflecting boundary conditions 
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Table 14. Numerical results of the c-r-4 scheme with non-reflecting boundary conditions 
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Figure 1 . — A surface element on the boundary 
S(V) of an arbitrary space-time volume V. 
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2(a) — The space-time mesh. 
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Figure 2. — The SEs and CEs. 
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Figure 3. — Definition of P , M , M + , and P + . 
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Figure 4.— A Burgers solution. 


REPORT DOCUMENTATION PAGE 

Form Approved 
OMB No. 0704-0188 

The public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 
data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this 
burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 
Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid OMB 
control number. 

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS. 

1. REPORT DATE (DD-MM-YYYY) 
01-08-2010 

2. REPORT TYPE 

Technical Memorandum 

3. DATES COVERED (From - To) 

4. TITLE AND SUBTITLE 

A New Approach for Constructing Highly Stable High Order CESE Schemes 

5a. CONTRACT NUMBER 

5b. GRANT NUMBER 



5c. PROGRAM ELEMENT NUMBER 

6. AUTHOR(S) 

Chang, Sin-Chung 

5d. PROJECT NUMBER 

5e. TASK NUMBER 

5f. WORK UNIT NUMBER 

WBS 599489.02.07.03.04.03.01 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
John H. Glenn Research Center at Lewis Field 
Cleveland, Ohio 44135-3191 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

E-17391 

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 

National Aeronautics and Space Administration 
Washington, DC 20546-0001 

10. SPONSORING/MONITOR'S 
ACRONYM(S) 

NASA 

11. SPONSORING/MONITORING 
REPORT NUMBER 

NASA/TM-20 10-2 16766 

12. DISTRIBUTION/AVAILABILITY STATEMENT 

Unclassified-Unlimited 
Subject Category: 64 

Available electronically at http://gltrs.grc.nasa.gov 

This publication is available from the NASA Center for AeroSpace Information, 443-757-5802 

13. SUPPLEMENTARY NOTES 

14. ABSTRACT 

A new approach is devised to construct high order CESE schemes which would avoid the common shortcomings of traditional high order 
schemes including: (a) susceptibility to computational instabilities; (b) computational inefficiency due to their local implicit nature (i.e., at 
each mesh points, need to solve a system of linear/nonlinear equations involving all the mesh variables associated with this mesh point); (c) 
use of large and elaborate stencils which complicates boundary treatments and also makes efficient parallel computing much harder; (d) 
difficulties in applications involving complex geometries; and (e) use of problem-specific techniques which are needed to overcome stability 
problems but often cause undesirable side effects. In fact it will be shown that, with the aid of a conceptual leap, one can build from a given 
2nd-order CESE scheme its 4th-, 6th-, 8th-,... order versions which have the same stencil and same stability conditions of the 2nd-order 
scheme, and also retain all other advantages of the latter scheme. A sketch of multidimensional extensions will also be provided. 

15. SUBJECT TERMS 

CESE method; High-order schemes 

16. SECURITY CLASSIFICATION OF: 

17. LIMITATION OF 
ABSTRACT 

uu 

18. NUMBER 
OF 

PAGES 

72 

19a. NAME OF RESPONSIBLE PERSON 

STI Help Desk (email:help@sti. nasa.gov) 

a. REPORT b. ABSTRACT c. THIS 

U U page 

u 

19b. TELEPHONE NUMBER (include area code) 

443-757-5802 


Standard Form 298 (Rev. 8-98) 
Prescribed by ANSI Std. Z39-18 

































